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ABSTRACT 

We study the relative velocity of inertial particles suspended in turbulent flows and discuss implica- 
tions of our results for dust particle collisions in protoplanetay disks. We carry out a 512 3 simulation 
of a weakly compressible turbulent flow, evolving 14 species of particles with friction timescale, r p , 
covering the entire scale range of the flow. The Stoke number, St, of the smallest particles is ~ 0.1, 
where St is the ratio of r p to the Kolmorgorov timescale, while the largest particles have t p ~ 54TL, 
where XL is the Lagrangian correlation timescale of the simulated flow. A comparison with the re- 
sults of the simulation shows that the theoretical model by Pan & Padoan (hereafter PP10) gives 
satisfactory predictions for the rms relative velocity between identical particles. The model reveals an 
insightful physical picture. It shows that the relative velocity of two same-size particles is determined 
by the particle memory of the flow velocity difference along their trajectories, and hence has a crucial 
dependence on the particle pair separation backward in time. Using the simulation data, we compute 
the collision kernel accounting for the effect of turbulent clustering. The collision kernel per unit cross 
section shows an abrupt rise as St increases toward ~ 1, which may be viewed as an activation process 
corresponding the rapid formation of caustics. At St ^ 1, the collision kernel first increases slightly, 
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reaches a maximum at r p ~ 2Tl, and finally decreases as r p for r p 3> 7l. Coagulation models 
of dust particle growth in protoplanetary disks should account for the new features of the collision 
kernel found in this study. Furthermore, we show that the probability distribution function (PDF) of 
the particle relative velocity is highly non-Gaussian, exhibiting extremely fat tails. Using the physical 
picture of PP10, we identify two sources of non-Gaussianity: the imprint of the turbulent flow inter- 
mittency and an intrinsic contribution from the particle dynamics. The shape of the PDF is fattest 
for particles with St ~ 1. The PDF tails for particles with r p ~ 1 — 2TL are well described by a 4/3 
stretched exponential function, consistent with a prediction based on the physical picture of PP10. 
The PDF approaches Gaussian only in the extreme large-particle limit with r p J> 54Tl. The high 
non-Gaussianity in the relative velocity is expected to have important implications for the growth and 
evolution of dust particles in protoplanetary disks, as it determines the fractions of collisions leading 
to sticking, bouncing or fragmentation. 



1. INTRODUCTION 

The dynamics of particles of finite inertia suspended 
in turbulent flows is a fundamental problem with ap- 
plications ranging from industrial processes (e.g. spray 
combustion engines) to geophysical flows (e.g., atmo- 
spheric clouds) . The interaction between turbulence and 
particles has been studied to understand rain initiation 
in warm terrestrial clouds (e.g., Pinsky & Khain 1997; 
Falkovich, Fouxon,& Stepanov 2002; Shaw 2003), cloud 
evolution in the atmospheres of planets, cool stars and 
brown dwarfs (e.g., Rossow 1978; Pruppacher & Klett 
1997; Freytag et al. 2010; Helling et al. 2011), collisions 
and growth of dust particles in protoplanetory disks (e.g., 
Dullemond & Dominik 2005; Zsom et al. 2010, 2011; 
Birnsticl et al. 2011) and in the interstellar medium (e.g., 
Ormel et al. 2009). 

The evolution of the particle size depends on the parti- 
cle collision rate which may be significantly enhanced by 
turbulent motions in the carrier flow, as illustrated by 
recent numerical and theoretical advances in this field 
(e.g., Wang et al. 2000; Zhou et al. 2001; Falkovich et 
al. 2002; Zaichik k Alipchenkov 2003, 2009; Zaichik et 



al. 2003, 2006; Wilkinson et al. 2005, 2006; Falkovich & 
Pumir 2007; Gustavsson & Mehlig 2011; Gustavsson et 
al. 2012). It is now understood that an accurate evalu- 
ation of the collision rates requires the modeling of two 
enhancement mechanisms: the preferential concentration 
or clustering of inertial particles (e.g., Maxey 1987 and 
Squires & Eaton 1991) and the turbulence-induced colli- 
sion velocity. In this work, we will focus on the statistics 
of turbulence-induced relative velocities, and briefly dis- 
cuss the contribution to the collision rate by turbulent 
clustering (see Pan et al. 2011 for a detailed discussion of 
turbulent clustering in the context of planetesimal forma- 
tion). We restrict our discussion to the collision velocity 
between same-size particles, usually referred to as the 
monodisperse case, and will address the general bidis- 
perse case (collisions between particles of different sizes) 
in a follow-up paper. 

The main motivation of our study is to improve the 
modeling of the evolution of dust particles in protoplan- 
etary disks, which sets the stage for the formation of 
planetesimals, the likely precursors to fully-fledged plan- 
ets. For example, the planetesimal formation model 
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by Johansen et al. (2007, 2009, 2011) requires particle 
growth up to decimeter to meter size, in order to achieve 
good frictional coupling to the disk rotation and hence 
the maximum clustering effect by the streaming insta- 
bility. Cuzzi et al. (2008, 2010) 'and Chambers (2010) 
proposed an alternative model of planetesimal formation 
based on the strong turbulent clustering of chondrulc- 
size particles. Other studies (e.g. Lee et al. 2010) focus 
on the possibility that small particles settle to the disk 
midplane, where gravitational instability can result in 
planetesimal formation (e.g. Goldreich & Ward 1973; 
Youdin 2011), despite the turbulence stirring caused by 
the Kelvin-Helmholtz instability induced by the vertical 
settling of the particles (e.g. Weidenschilling 1980; Chi- 
ang 2008). 

The evolution of the size distribution of dust particles 
is controlled by collisions. Small particles tend to stick 
together when colliding, and thus their size grows by co- 
agulation. As the size increases, the particles become 
less sticky (Blum & Wurm 2010), and, depending on the 
collision velocity, the collisions may result in bouncing or 
fragmentation. A detailed summary of experimental re- 
sults for the dependence of the collision outcome on the 
particle properties (such as the particle size and poros- 
ity) and on the collision velocity can be found in Guttler 
et al. (2010). The coagulation, bouncing and fragmenta- 
tion processes may lead to a quasi-equilibrium distribu- 
tion of particle sizes (e.g., Birnstiel et al. 2011; Zsom et 
al. 2010, 2011). Due to the dependence of the collision 
outcome on the collision velocity, an accurate evaluation 
of the turbulence-induced relative velocity is important 
for modeling the size distribution of dust particles. 

Saffman and Turner (1956) studied the relative veloc- 
ity in the limit of small particles with the particle friction 
or stopping time, r p , much smaller than the Kolmorgorov 
timescale, t v , of the turbulent flow. This limit, known as 
the Saffman- Turner limit, is usually expressed as St <C 1, 
where the Stokes number is defined as St = t p /t v . 
Saffman and Turner (1956) predicted that, at a given 
particle distance, r, the relative velocity of identical par- 
ticles is independent of St, and, at a given St, it scales 
linearly with r for small r. In the opposite limit of large 
particles with r p larger than the largest timescale of the 
turbulent flow, Abrahamson (1975) showed that the rel- 
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ative velocity scales with the friction time as r p ' . A 
variety of models have been developed to bridge the two 
limits and predict the relative velocity for particles of any 
size, i.e., with r p covering the entire scale range of the 
carrier flow (Williams & Crane 1983, Yuu 1984, Kruis & 
Kusters 1997, & Alipchenkov 2003, Zaichik et al. 2006, 
Ayala et al. 2008). Among these models, the formulation 
of Zaichik and collaborators is particularly impressive, 
as it examines both turbulent clustering and turbulence- 
induced relative velocity simultaneously. The model pre- 
diction for the relative velocity agrees well with simula- 
tion results at low resolutions. However, the model lacks 
a transparent physical picture. 

Pan & Padoan (2010) developed a new model for the 
relative velocity of inertial particles of any size that pro- 
vides an insightful physical picture of the problem. Their 
formulation illustrates that the relative velocity of iden- 
tical particles is determined by the memory of the flow 
velocity difference along their trajectories in the past. 



The model also shows that the separation of inertial par- 
ticle pairs backward in time plays an important role in 
the relative velocity of nearby particles. The model can 
correctly reproduce the scaling behaviors of the relative 
speed in the extreme limits of small and large particles 
(see above), and its prediction was found to successfully 
match the simulation data of Wang et al. (2000). In this 
paper, we will further test the model against numerical 
simulations of considerably higher resolution. 

Falkovich et al. (2002) discovered an interesting effect, 
named the sling effect, which provides an important con- 
tribution to the collision rate. The basic physical picture 
of this effect is that inertial particles may be shot out of 
fluid streamlines with high curvature, causing their tra- 
jectories to cross with those of other particles (Falkovich 
& Pumir 2007). In particular, in flow regions with large 
negative velocity gradients, fast particles can catch up 
with the slower ones from behind. For small particles 
with St -C 1, the sling events correspond to high-order 
statistics of the flow velocity gradient, and the effect is 
not reflected in the prediction of Saffman and Turner 
(1956). The trajectory crossing causes the particle ve- 
locity to be multi- valued at a given point, leading to the 
formation of caustics in the momentum-position phase 
space of the particles (Wilkinson et al. 2006). The slings 
or caustics increase the local particle density, and give 
rise to a finite and significant collision velocity for par- 
ticles even at an infinitesimal distance (Falkovich et al. 

2002, Wilkinson et al. 2006). Both effects tend to en- 
hance the particle collision rate. It has been shown that, 
as St approaches 1, the effect of slings or caustics causes 
a rapid rise in the collision rate, which has been pro- 
posed to be responsible for the initiation of rain shower 
in terrestrial clouds (Wilkinson et al. 2006). Applying 
this effect to dust particle collisions in protoplanetary 
disks, one may expect that the collision rate greatly ac- 
celerates as the particle grows past sub-mm to mm size, 
corresponding to St ~ 1 for typical protoplanetary tur- 
bulence conditions. 

The recent developments mentioned above have not 
been considered in coagulation models for dust particles 
in circumstcllar disks. We find that, in the astrophysical 
literature, the general formulation of the collision kernel 
of dust coagulation models is inaccurate. In particular, 
the formulation significantly underestimates the collision 
rate for particles with St ~ 1, because it ignores the ef- 
fect of turbulent clustering. Furthermore, the dust co- 
agulation models usually adopt collision velocities from 
the work of Volk et al. (1980) and its later extensions 
(e.g., Markiewicz, Mizuno & Volk 1991, Cuzzi and Hogan 

2003, and Ormel & Cuzzi 2007), which has a number of 
severe limitations. First, Pan & Padoan (2010) pointed 
out a weakness in the physical picture of these models. 
Roughly speaking, these models assume the velocities of 
two particles induced by turbulent eddies with turnover 
time significantly smaller (larger) than r p are indepen- 
dent (correlated). As shown by Pan & Padoan (2010), 
whether the particle velocities contributed by turbulent 
eddies of a given size are correlated or not also depends 
on how the eddy size compares to the separation of the 
particles at the time the eddies were encountered, so the 
eddy turnover time is not the only factor that determines 
the degree of correlation. The role of the particle separa- 
tion relative to the eddy size is not captured by the ap- 
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proach of Volk et al. Second, the models do not account 
for the dependence of the relative velocity on the parti- 
cle distance, r, and, as a result, they cannot reproduce 
the Saffman- Turner regime for small particles. We will 
show that the relative velocity does have a r-dependence 
for St ^ 6 particles, which, together with the scale de- 
pendence of turbulent clustering, is crucial for the eval- 
uation for the collision kernel of these particles. Finally, 
the model of Volk et al. is found to overestimate, by a 
factor of 2, the relative velocity of particles with r p of 
the order of the large eddy turnover time. 

In this paper, we use numerical simulations to study 
inertial particle dynamics in turbulent flows. With the 
simulation data, we test the model of Pan & Padoan 
(2010) and compute the collision kernel as a function of 
St. We also explore the probability distribution func- 
tion (PDF) of the relative velocity, which is expected to 
have important effects on the particle size distribution. 
As mentioned earlier, the outcome of dust particle colli- 
sions depends on the collision velocity, and thus, due to 
the stochastic nature of turbulence-induced relative ve- 
locity, collisions between particles of the same properties 
may lead to different outcomes. The PDF of the collision 
velocity is needed to estimate the fractions of collisions 
leading to sticking, bouncing or fragmentation. The PDF 
of the particle relative velocity has been shown to be 
highly non-Gaussian by various numerical, experimental 
and theoretical studies (e.g., Sundaram and Collins 1997, 
Wang et al. 2000, Gustavsson et al. 2008, Bee et al. 2009, 
de Jong et al. 2010, Gustavsson & Mehlig 2011, Hubbard 
2012). This highly non-Gaussian nature of the relative 
velocity should be incorporated into coagulation models 
for the size distribution of dust particles in protoplanc- 
tary disks. 

In this work, we will investigate the particle dynam- 
ics only in statistically homogeneous and isotropic tur- 
bulence. This is clearly an idealized situation, consid- 
ering various complexities in protoplanetary disks. For 
example, the disk rotation induces large-scale anisotropy, 
which may have significant effects on the prediction for 
particles with friction time close to the rotation period. 
Nevertheless, the idealized problem is a very useful tool 
to understand the fundamental physics, and it provides 
the first step toward accurately modeling the dust par- 
ticle collision velocity in protoplanetary turbulence. We 
also neglect the vertical settling and radial drift. These 
processes do not directly affect the relative velocity be- 
tween identical particles, although they probably provide 
important contributions for particles of different sizes 
that we address in a follow-up work. 

The paper is organized as follows. In §2, we present 
a simple model for the rms velocity of a single particle, 
which provides an illustration for our formulation of the 
particle relative velocity. In §3, we introduce the model 
of Pan & Padoan (2010) for the relative velocity of nearby 
particles, review the physical picture of the model, and 
discuss the qualitative behavior of the model predictions. 
Our simulation setup and the statistical properties of the 
simulated turbulent flow are described in §4. §5 presents 
simulation results for the particle rms velocity. In §6, we 
test the model of Pan & Padoan (2010) with respect to 
the rms relative velocity of particles found in the simula- 
tion, compute the collision kernel accounting for the ef- 
fects of both turbulent clustering and turbulence-induced 



collision velocities, and present in detail the probability 
distribution of the relative velocity as a function of the 
particle inertia. The conclusions of our study are sum- 
marized in §7. 

2. THE VELOCITY OF INERTIAL PARTICLES 

The dynamics of inertial particles depends crucially on 
its friction or stopping timescale, r p . To evaluate of the 
friction timescale, we first need to compare the particle 
size, a p , with the mean free path of the gas particles 
in the carrying flow. If the particle size is larger than 
the mean free path, the friction timescale is given by the 

Stokes law r p = § (^j, where p d (~ 1 g cm" 3 ) 

is the density of the dust material, p is the gas density, 
and v is the kinematic viscosity of the flow. On the other 
hand, if a p is smaller than the gas mean free path, the 

particle is in the Epstein regime and r p = f^fj (^)' 
where C s is the sound speed in the flow. For example, 
for a typical value of the gas density in protoplanetary 
discs, p ~ 10~ 9 g cm~ 3 at a radius of 1 AU, the mean 
free path of the gas particles is about 1 cm, and thus 
particles with a p larger (smaller) than 1 cm are in the 
Stokes (Epstein) regime. 

The velocity, v(t), of an inertial particle suspended in 
a turbulent velocity field, u(x,t), obeys the equation 

dv u (X(t), t) — v 

j. i (1) 

at t p 

where X(t) is the position of the particle at time t, and 
u (X(t), t) corresponds to the flow velocity "seen" by the 
particle. Eq. ([1]) has a formal solution, 

v(t) = - f u (X(t),t) exp (-—) dr, (2) 

T P J t V T P / 

where it is assumed that t — 1 >• t p and the particle has 
already lost the memory of its initial velocity at to. The 
formal solution indicates that the velocity of an inertial 
particle is determined by the memory of the flow velocity 
along its trajectory within a timescale of ~ t p in the past. 

Although the aim of the present work is the collision ve- 
locity of inertial particle pairs, we start with a discussion 
of the single-particle (or "1-particle") velocity induced 
by turbulent motions. We provide a simple model for 
the 1-particle rms velocity as a function of the particle 
friction time. The derivation of this model helps to illus- 
trate our formulation for the collision velocity between 
two nearby particles. 

The 1-particle rms velocity can be calculated using the 
formal solution of eq. ([2]) . We assume the turbulent flow 
is statistically stationary, and the particle statistics even- 
tually relax to a steady state. We consider a time when 
the steady state is already reached and, without loss of 
generality, we denote this time as time 0. Using eq. ^ 
at t = 0, we have, 

f° dr f° dr' fr\ (t>\ 

{ViVj) = — / — B Tij [t, t ) exp I — 1 exp I — 1 , 

J - oo T p J - oo T p \ T P / \ T P J 

(3) 

where B Tij (r,T') = (ttj(X(r), r)uj (X(t'), t')) is the 
temporal correlation tensor of the flow velocity along the 
trajectory, X(t), of the inertial particle. The subscript 
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"T" stands for "trajectory" . We have also changed the 
lower integration limit (to) in eq. ([2]) to — oo, based on 
the assumption that the particle dynamics has been fully 
relaxed at time (i.e., to <C — r p ). 

With statistical stationarity and isotropy, the trajec- 
tory correlation tensor can be written as Bnj(r,T') = 
u' 2 $i(t' — T~)5ij, where u' is the ID rms velocity of the 
turbulent flow and the correlation coefficient $i is a func- 
tion of the time lag only. The subscript "1" is used to 
indicate that the correlation is along the trajectory of 
one particle. The correlation coefficient, $i, is unknown, 
and a common assumption is to approximate it with the 
Lagrangian correlation function, $l, of tracer particles 
(or fluid elements), which has been extensively studied 
in the literature. The assumption is likely valid for small 
particles, but cannot be justified for large particles on 
a theoretical basis. We will validate the assumption a 
posteriori using simulation results. 

The simplest choice for the form of $l is an ex- 
ponential function, $l(At) = exp(— |At|/2l), where 
At = t' — t is the time lag and Tl the Lagrangian corre- 
lation timescale. Setting B-nj — u' 2 exp(— |r' — t|/Tl)% 



in eq. ([3]), we have 
particle velocity, v', is given by, 



~Sij, where the ID rms 



T l + t p 



1/2 



(4) 



This result shows that the particle rms velocity ap- 
proaches the flow velocity for r p <C Tl and decreases 
as (Tl/tp) 1 / 2 for r p > T L (e.g., Abrahamson 1975). In 
the large particle limit, t p ;g> Tl, the action of even the 
largest turbulent eddies on the particle would appear to 
be random kicks when viewed on a timescale of t p . In 
that case, eq. ([1]) is essentially a Langevin equation, and 
the particle motions are similar to Brownian motions. 

— 1/2 

The t p scaling corresponds to an "equilibrium" be- 
tween the velocity of these large particles and the turbu- 
lent motions of the flow. 

Numerical simulations have shown that the Lagrangian 
correlation function, $l(At), is better fit by a bi- 
exponential form (e.g., Sawford 1991). A single- 
exponential form does not reflect the smooth part of the 
correlation function for time lag At smaller than the 
Taylor micro timescale, Tp- The Taylor timescale is de- 

1/2 

fined as (2u' 2 /a 2 ) , where a is the rms acceleration of 
the turbulent velocity field. The bi-exponential form for 
$l(At) is, 



$l(At) 



1 



2-s/l - 2z 2 



(1 + v 7 ! - 2z 2 ) x 



exp 



exp 



2|At| 



(l + Vl-2z 2 )T L 
2|At| 



- (l - \/l-2z 2 ): 



l-yT^2F T L 



(5) 



where the parameter z is defined as the ratio, z = 
tt/Tl- From the above equation, it is easy to show that 
Tl = / $l(At)c?At, and the bi-exponential function is 



smooth, ~ 1 — (At/tt) 2 , at At <C Tt- In the limit 
z — > 0, eq. ([5]) reduces to the single exponential with a 
timescale of Tl. 

Adopting the bi-exponential form, eq. ([5]), for the tra- 
jectory correlation coefficient, $ 1; we find that the one- 
particle rms velocity is given by, 



n + z 2 /2 
n + n 2 + z 2 /2 



1/2 



(6) 



where 57 is defined as 57 = t p /Tl. In the limits 57 <C 1 
and f! > 1, eq. ((6]) has the same behavior as eq. (j4j 
from the single exponential correlation. In fact, the two 
predictions, eqs. (U]) and ©, are close to each other 
at all values of 57, differing only by a few percent at 
57 ~ 1. This suggests that, for a given correlation 
timescale, Tl(= / <&l(At)g?At), the integral in eq. ((3]) 
is insensitive to the exact function form of the correla- 
tion, $i(At). We will measure the parameters, z and Tl, 
using Lagrangian tracer particles in our simulated turbu- 
lent flow, and test the predictions, eqs. (fj| and ([5]), for 
the 1-particle rms velocity against the simulation data. 

3. TURBULENCE-INDUCED RELATIVE VELOCITY OF 
INERTIAL PARTICLES 

We briefly review the 2-point Eulerian statistics of 
the velocity field in fully-developed turbulence, which 
is crucial to understand the relative velocity of two in- 
ertial particles. We consider the 2nd-order structure 
tensor, defined as Sij(£) = (AuiAuj) where Aui — 
Ui(x + £, t) — Ui(x,t) is the velocity increment across 
a separation £ and (• • •) denotes the ensemble average. 
The statistics of Au; is independent of x and t from the 
assumption of homogeneity and stationarity. With sta- 
tistical isotropy, the velocity structure tensor takes the 
form (e.g., Monin and Yaglom 1975), 



S ij (£) = S Dn (£)8 ij + [Sn{i) 



Snn(£)] —pT 



(7) 



where the longitudinal and transverse structure func- 
tions, Su and iSnn, are functions of the amplitude, £, of 
the separation, £, but not of its direction, £/£. From 
eq. ©, we see Su = Sij(£)£i£j/£ 2 = (Au 2 ), where 
Ait r = Aui£i/£ is the radial component of the flow veloc- 
ity increment. Similarly, the transverse structure func- 
tion S nn can be written as S nn = ((Au t ) 2 ) with Au t 
being one of the two components of Am on the tangen- 
tial/transverse plane perpendicular to £. The statisti- 
cal isotropy indicates that the probability distribution 
of Aut is invariant under any rotation about the direc- 
tion £/£. In incompressible turbulence, which is approxi- 
mately the case for gas flows in protoplanetary disks, we 
have the relation S nn = Su + ^idS\\/d£, which can be de- 
rived from the incompressibility condition: djSij(£) = 
(Monin and Yaglom 1975). 

The structure functions exhibit different scaling behav- 
iors in different scale ranges. There are three subranges 
divided by two length scales, the Kolmorgorov length 
scale, 77, and the integral length scale L. The Kolmor- 
gorov scale, T], is defined as r\ — (i^/e) 1 / 4 , where v and 
e are, respectively, the kinematic viscosity and the av- 
erage energy dissipation rate in the turbulent flow. It 
essentially corresponds to the size of the smallest eddies. 
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Scales below r\ are called the viscous or dissipation range, 
where the velocity field is laminar and differentiable due 
to the smoothing effect of the viscosity. In the dissi- 
pation range (£ rf), the velocity difference scales lin- 
early with £, and the longitudinal structure function is 
Sn = yf^ 2 . S nn is twice larger, i.e., S nn = j^£ 2 , as re- 
quired by the incompressibility constraint. In the inertial 
range, r\ <, I < L, Sn follows the Kolmorgorov scaling, 

5*11 = Ck^) 2 / 3 , where the coefficient Ck is known as the 
Kolmorgorov constant. The typical value of Ck is ~ 2. 
The incompressibility condition gives S lln = 45n/3 in the 
inertial range. The integral scale, L, is essentially the 
correlation length scale of the velocity field. At I ^> L, 
the velocity field is statistically uncorrelated, and both 
S\i(£) and S nn (£) are constant and equal to 2u' 2 with u' 
the ID rms velocity of the flow. 

To bridge the scalings of 5*11 in the three scale ranges, 
we adopt a connecting formula (Zaichik et al. 2006), 



1 — exp 



4/3 



(15C K ) 3 / 4 

(i/v) 4 



-1 1/6 



(£/ V y + (2u' 2 /C K u 2 f 



(8) 



where u v is the Kolmorgorov velocity scale defined as 
(Ve) 1 / 4 . With eq. ([8]) for Sn, we can obtain S nn using the 
incompressibility condition (see above). Alternatively, 
one may adopt a separate connecting formula for S nn . 
In §4.3, we will measure the structure functions in our 
simulated flow, and fit them with connecting formulas. 

The goal of this work is to understand and predict the 
velocity at which two nearby inertial particles collide. 
The collision speed of two particles is essentially their 
relative velocity as they cross, i.e., as their distance, r, 
becomes equal to the sum of the particle radii (Saffman 
and Turner 1956). For applications to dust particles in 
protoplanetary disks, the particle size is typically much 
smaller than the Kolmorgorov length scale, rj. In the 
following, we therefore explore the relative speed of par- 
ticles at small distances, r <S r\. 

We label two particles coming together with super- 
scripts (1) and (2). For example, we denote their po- 
sitions as X^(t) and X^'{t), and their velocities as 
v^\t) and v^(t) (see Fig. 1 for illustration). When 
the superscripts (1) and (2) are not used, the discus- 
sion is general and not referring to a specific particle. 
At a given time t, we consider the relative velocity w = 
v( 2 \t) — «W(t), of particle pairs at a given separation, r, 

which corresponds to a constraint X^ 2 '(t) — X^'(t) = r 
for the particle positions. We first present a theoreti- 
cal model for the second-order moment of w, and then 
use numerical simulations to explore its full statistics, in- 
cluding the probability distribution function (PDF), as 
a function of r p . 

Similar to the structure tensor of the flow velocity, we 
characterize the second-order statistics of the particle rel- 
ative velocity w by a structure tensor, 



(9) 



which was referred to as the particle velocity structure 
tensor by Pan and Padoan (2010). The tensor S P ij can 



be evaluated using the formal solution, eq. ©, of the par- 
ticle velocity and the statistical properties of the carrier 
flow (see Pan and Padoan 2010 and §3.2). 

Once the particle dynamics is fully relaxed, the particle 
velocity is expected to possess the same statistical sym- 
metries as the flow, including stationarity, homogeneity 
and isotropy. With these symmetries, the particle struc- 
ture tensor S P ij can be written in a similar form as the 
structure tensor of the flow (eq. [7]), 

S pij (r) = K 2 )<5 y - + (K 2 ) - K 2 )) ^, (10) 



where 



and 



are the variances of the ra- 



dial/longitudinal component, w T (= WiTi/r), and a tan- 
gential/transverse component, wt, of the relative veloc- 
ity, respectively. For particle collisions, we are only inter- 
ested in S P ij at r ^ r\. Under the assumption of isotropy, 
the tangential component, u>t, is expected to be statis- 
tically invariant for any rotations about the axis r. We 
can thus measure the statistics of the tangential relative 
velocity by projecting w into an arbitrary direction on 
the plane perpendicular to r. 

The radial relative speed, u> r , is of particular inter- 
est for the estimate of the particle collision frequency 
(Wang et al. 2000; see §6.2), and its variance is related 
to the particle structure tensor as (w 2 ) = Spijrirj/r 2 . 
The total collisional energy of a particle pair plays a cru- 
cial role in the outcome of the collision (e.g., Blum and 
Wurm 2008). The average collisional energy can be com- 
puted from the 3D variance, (w 2 ), of the relative velocity, 
which is the contraction of the particle structure tensor, 
i.e., Spu = (w 2 ) + 2(w 2 ). In the rest of this section, 
we consider theoretical models for the variances of the 
relative velocity. As mentioned earlier, we focus on the 
monodisperse case with equal-size particles. 

3.1. The Limits of Small and Large Particles 

We first consider small particles in the Saffman- Turner 
limit (hereafter the S-T limit). In this limit, the fric- 
tion timescale, r p , is much smaller than the Kolmorgorov 
timescale, t^, of the carried flow, which is defined as 
T n = {v/e) 1 / 2 . The Kolmor gorov timescale is the small- 
est timescale in a turbulent flow, corresponding to the 
turnover time of the smallest eddies. Therefore, the ve- 
locity of particles with t p <ti can be approximated by 
a Taylor expansion of eq. ([1}, v (t) — u(X,t)+T p a(X,t), 
where a = Du/Dt is the acceleration of the local fluid 
element. Applying the approximation to both particles 
(1) and (2), we have w = (u^ - u^) + (a^ - a^) r p , 

where m^ 2 ) (= u{X^ 2 \t)) and a^ 2 ) (= a{X^ 2 \t)) 
are the flow velocity and acceleration at the positions of 
particles (1) and (2), respectively. Saffman and Turner 
(1956) assumed that the correlation coefficient of the flow 
accelerations, and a^ 2 \ across a small distance, r, is 
unity, which is equivalent to assuming ~ a^ 2 \ The 
acceleration terms then cancel out for identical particles, 
and the particle structure tensor, S P ij, is simply equal 
to the flow structure tensor, Sij, defined in eq. ([7]). Us- 
ing the flow structure functions Sn and S^n at r <ti rj in 
incompressible turbulence, we have the Saffman- Turner 
formula, 



/ 2\ 1 6 2 / 2\ 2 e 2 



(11) 
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for identical particles with St <C 1. The equation shows 
that in the S-T limit the relative speed is caused by the 
flow velocity difference across the particle separation. 
The effect is usually referred to as the shear contribu- 
tior£J From eq. flTTl) . the 3D variance of the relative 
velocity is given by (w 2 (r)) = -^r 2 . 

The S-T formula predicts that the tangential variance 
of the relative velocity, (w 2 ), is twice larger than that 
in the radial direction. Eq. (|f f [) also indicates a con- 
stant relative speed at a given separation, r, and a lin- 
ear scaling with r at a given St <C 1. The accuracy of 
the Saffman- Turner formula has been questioned, which 
neglects the effect of slings and caustic formation (e.g., 
Falkovich et al. 2002, Wilkinson et al. 2006). We will test 
the S-T prediction against our simulation data. In the 
S-T limit, the particle memory is short and the relative 
speed is determined largely by the local flow velocity at 
small scales. The memory effect becomes important for 
larger particles with r p > t v (see §3.2). 

We next consider the other extreme limit, i.e., large 
particles with r p much larger than the Lagrangian corre- 
lation time, Tl, of the flow. As discussed in §2, the mo- 
tions of these particles are similar to Brownian motions, 
and the velocities of any two particles are statistically in- 
dependent. This is because the velocity of a large particle 
has a significant contribution from its memory of the flow 
velocity long time ago, and the flow velocities "seen" by 
the two particles at that time were uncorrelated because 
the particle separation was likely larger than the flow in- 
tegral length scale, L. With the independence of v^ 1 ' 
and v^ 2 \ the particle structure tensor defined in eq. ^ 



can be written as S P ij = 



y(2)l 



Sij, where 



u'W and v'( 2 ) are the (ID) rms velocities of particles (1) 
and (2), respectively. As shown in §2, for particles with 

1/2 

r p ^> Tl, the rms velocity is given by ~ v! (Tl/t p ) 
We therefore have (e.g., Abrahamson 1975), 



= 2u' 2 ^ 



(12) 



for identical particles with r p ^> Tl- The equation sug- 
gests that the particle relative speed decreases with the 
Stokes number as St^ 1 / 2 . The physical picture for the 
relative velocity in the large particle limit is clear, and 
the result, eq. (fl2j) . is robust. 

In between the two extreme limits are particles in 
the inertial range, i.e., particles with friction timescale 
t v Ss t p <, Tl, corresponding to inertial-range scales in 
the turbulent flow. The physics for the relative velocity 
of inertial-range particles is more complicated. For ex- 
ample, in the two extreme limits the velocity correlation 
of two nearby particles is relatively easy to estimate: The 
velocities are either highly correlated (for small particles) 
or essentially independent (for large particles). On the 
other hand, for inertial-range particles, the velocity cor- 
relation is at an intermediate level. We will show that 
the key physics for the relative velocity of inertial-range 

2 The term "shear contribution" is as opposed to the "acceler- 
ation contribution" from the acceleration terms mentioned above, 
which do not vanish for particles of different sizes. The acceleration 
contribution in the bidisperse case will be discussed in a separate 
paper. 



particles is the particles' memory of the flow velocity dif- 
ference in the past and the separation of the particle pair 
backward in time. 

As mentioned in the Introduction, a variety of mod- 
els for the particle relative velocity covering the whole 
range of particle sizes have been developed (e.g., Volk et 
al. 1980, Ormel & Cuzzi 2007, Zaichik & Alipchenkov 
2003, Zaichik et al. 2006, and Pan & Padoan 2010). The 
models listed here all predict a St 1 / 2 scaling for inertial- 
range particles in turbulent flows with an extended in- 
ertial range. The models of Zaichik and collaborators 
and Pan and Padoan (2010) can reproduce both the S-T 
limit (eq. JTTJ)) and the large particle limit (eq. ([T21Q. We 
will focus on the model of Pan and Padoan (2010), which 
provides a clearer physical picture than that of Zaichik 
et al. The physical differences between various models 
have been summarized in Pan and Padoan (2010). 

3.2. The Model of Pan and Padoan (2010) 

We first review the formulation and the physical pic- 
ture of the model by Pan & Padoan (2010; hereafter 
PP10) for the relative velocity of identical particles. The 
main idea of the PP10 model was to evaluate the par- 
ticle velocity structure tensor, S V ij{r), using the formal 
solution (eq. @) for the particle velocity. Applying eq. 
^ to the velocities of particles (1) and (2) at t = 0, we 
have, 



,(!) 



1 



a 



(1 



(r) 



dr, 



where «,( 1:2 )(t) (e 



exp 

\ t p/ 

(13) 

u(X^ 1,2 ' (t), t)) are the flow veloc- 
ities "seen" by the two particles at time r. Again we 
have changed the lower integration limit in the formal 
solution, eq. ©, to — oo. 

Inserting eq. (|13[) into the definition, eq 
is straightforward to find that, 



of Spij , it 



S p ij(r) 



dr 



S T ij(r;T,T)x 



exp 



— ) exp ( — 



(14) 



where Snj, named the trajectory structure tensor by 
PP10, is defined as, 



S Tij (r,T,T') = 



,(2) 



(t)- U W(t) 



,( 2 )/ 



V)-«5V)])- 

(15) 

Clearly, this tensor corresponds to the correlation of the 
flow velocity differences along the trajectories of the two 
particles at two times r and r' . The trajectory structure 
tensor, Srij, depends on the separation, r, through the 

constraint that X (2) (0)-X (1) (0) = r. Note that eq. JUJ 
is in close analogy with eq. ^ for the 1-particle velocity. 
Here the trajectory structure tensor, Sry, replaces the 
trajectory correlation tensor, B^ij, in eq. ([3]). 

The physical meaning of eqs. (fT3"|) and (fT4")) is clear: 
the relative velocity of two identical inertial particles is 
controlled by the particles' memory of the flow velocity 
difference within a friction timescale, ~ r p , in the past. 
The physical picture is illustrated in Fig. 1. The tra- 
jectory structure tensor, Stij, is unknown, and we next 
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model it using the approach of PP10. 

Since the flow velocity difference scales with the dis- 
tance, the trajectory structure tensor, S^ij, has an in- 
direct dependence on the particle separation at r and 
t'. Wc denote the particle separation at r as d(r) 

(= X (2) (t) - X (1) (t)). The separation vector d is 
stochastic because of the random dispersion of the parti- 
cle pair by turbulent motions. Stij also has a dependence 
on the time lag (V — r). This dependence is associated 
with the temporal correlation of turbulent structures or 
eddies encountered by the two particles between r and 
t', and the correlation time is essentially the turnover 
time of these eddies. To estimate Sry, we consider the 
(indirect) spatial dependence on the particle separation 
and the temporal dependence on the time lag separately. 

We use a typical particle separation R(t, t') between 
t and t' to model the spatial dependence. Clearly, 
R(t,t') is a random vector because of the stochastic- 
ity of d(r) and d(r'). We approximate the dependence 
on the separation by the Eulerian structure tensor of the 
flow velocity, Sij(R), defined in eq. ([7]). We denote as 
$2 (t' — r, R) the temporal correlation of the flow struc- 
ture at the scale R. The correlation, $2, is expected to 
be an even function of the time lag and is normalized 
to unity, $2(0, R) = 1, at zero time lag. To distinguish 
from the temporal correlation, $i(r' — r), along the tra- 
jectory of a single particle (see §2), we have used a sub- 
script "2" here for the two-particle case. The function 
form of $2 will be specified in §3.2.1. The trajectory 
structure tensor is then modeled as the product of the 
two dependences (PP10), 

S Tij (r;r,r') ~ (S^R^r' - r, R)) R (16) 

where (•••)« denotes the average over the statistics of the 
random vector, R. The average is over the probability 
distributions of both the amplitude, i?, and the direction 
of R. Eq. (TT6"f implicitly assumes the statistical inde- 
pendence of the velocity difference, Am, seen by the two 
particles from their separation, R. Rigorously, the ampli- 
tudes of Am and R may have a correlation. If the particle 
pair encounters an eddy with a larger velocity, the par- 
ticle separation tends to be larger. For example, if R is 

in the inertial range of the flow, Au ~ e H /3 i? 1/3 from the 
refined similarity hypothesis (Kolmorgorov 1962), where 
en is the average dissipation rate over the scale R as seen 
by the particle pair. A positive correlation is expected 
between eR and R. Eq. ([To]) neglects this correlation and 
may underestimate St^ and hence the particle relative 
velocity. 

The $2 term in eq. (| 10[) does not depend on the direc- 
tion of R, so one can first take the angular average of 
Sij(R) and then average the entire term over the PDF 
of the amplitude, R. The latter cannot be exactly per- 
formed because the PDF of R is unknown. With some 
simple estimates, PP10 argued that simply using the rms 
of R to evaluate Sxy instead of averaging over the PDF 
of R only causes a small difference (~ 10%) in the model 
prediction. Following PP10, we ignore the PDF of R and 
insert the rms of R to evaluate 5t»j. For the simplicity 
of notations, we use R to denote the rms of the particle 
distance in between r and r' in the rest of the paper. A 
similar notation is adopted for d(r) and d(r'), which will 
denote the rms particle separations at r and t', respec- 



tively. We approximate R by the geometric average of 
d(r) and d(r'), 

R(T,r') = [d(T)d(r')} 1/2 . (17) 

The rms separation d(r) as a function of time r will be 
discussed in §3.2.3. 

With the above assumptions, the trajectory structure 
tensor is modeled as, 

S Tij (r;r,r') ~ (S tJ (R)) ^{t' - r, R). (18) 

The angular average over the direction of R will be car- 
ried out in §3.2.2. In eq. ([15]). the dependence of Sr-y on 
r is through the dependence of the particle separations, 
d(r), d(r') and R, on r. We will refer to r as the "ini- 
tial" separation, although it actually corresponds to the 
current or final separation of the two particles. Clearly, 
our formulation indicates that the separation of particle 
pairs backward in time is crucial for the prediction of the 
particle relative speed. 

Inserting eq. (fT8)l into eq. (|T4l) gives the PP10 model 
for the particle structure tensor, 




$ 2 (T'-T,i?)expf-^ exp(^) . (19) 

We will numerically compute this double integral after 
evaluating or modeling the angular average, the temporal 
correlation function and the particle separation backward 
in time. 

A simplification of the PP10 model is to set R to one 
of two distances, d(r) or d(r'), instead of their geometric 
average. We find that replacing R in eq. (fT9|) by either 
d(r) or d(r') leads to equivalent model prediction for 
the particle relative speed. This is because the temporal 
correlation $2 in eq. (|T9")) is an even function of Ar(= 
t' — r), and the product of the two exponential cutoffs 
are invariant under the exchange of r and r 1 . If one sets 
R = d(r) in eq. (fT9|) . the integral over r' can be isolated, 
yielding, 

S plJ = -[ (Sij(d(j)))^F{T)exp(—)dT, (20) 
t p J -00 \ t pJ 

where the angular average is over the direction of d(r) 
and the function F(t) is defined as, 

F(r) = - [ $ 2 (r' - r,d(r)) exp (-) dr 1 . (21) 

The factor F[t) may be roughly viewed as a response 
function of the particle pair to turbulent eddies at the 
scale d(r). Although not indicated explicitly, the factor 
F(t) also depends on the initial particle separation, r, 
through its dependence on d(r). We will refer to eqs. 
(|2"D|) and (l2~Tj) as the simplified model. The advantage of 
the simplified model is that the response function F(t) 
can be integrated analytically using assumed function 
forms of $2 in §3.2.1, and one only needs to numerically 
solve a single integral in eq. (f2"0)l . On the other hand, for 
the original PP10 model, one must numerically evaluate 
the double integral in eq. (|T9| . 
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Particle 2 
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Fig. 1. — Schematic figure illustrating the physical picture of the PP10 model for the relative velocity of inertial particles of equal size. 
At time 0, the separation of the two inertial particles, (1) and (2), is r. The velocity, v(0), of each particle at t = 0, is determined by 
its memory of the flow velocity, it(r), along the particle trajectory in the past. The relative velocity of the two particles depends on 
the flow velocity difference, u^ 2 \t) — u^'(t), they "saw" within about a friction timescale t p in the past, i.e., — t p < r < 0. The flow 
velocity difference at a given time r is related to the particle separation, d(r). The separation of the two particles satisfies the "initial" 
constraint d(0) = r and increases backward in time. Due to particle inertia, a roughly ballistic separation is expected within a particle 
friction timescale. The trajectories plot here reflect a more-or-less linear separation of the two particles. The temporal correlation of the 
flow velocity differences the two particles "saw" at different times, say r and r', is another important effect The correlation timescale is 
associated with the turnover time for turbulent eddies encountered by the two particles between r and r'. 



3.2.1. The temporal correlation $2 

To estimate the temporal correlation in the trajectory 
structure tensor, Spy, we first consider a special case 
where the typical particle separation, R, is much larger 
than the integral length scale, L, of the flow. In this 
case, the flow velocities, u^> and it' 2 ', "seen" by the 
two particles are independent, and Sj^j can be writ- 

ten as (u{ 1) (t)uJ 1) (t / )> + (uf ] ( T )uf 3 (r')) (see eq. ©). 
Both terms correspond to the trajectory correlation ten- 
sor B^ij defined below eq. ([3]) in §2, and for identical 
particles the two terms are equal. Therefore, for R 3> L, 
the correlation function $2 (At, R) is the same as the 
temporal correlation coefficient, $ 1 (Ar), along the tra- 
jectory of one particle. 

In §2, we approximated $1 by the Lagrangian corre- 
lation function, $l- Using the approximation again, we 
have $ 2 (At,£) = $i(Ar) ~ $ L (Ar) for £ > L. Two 
function forms, single- and bi-exponential, were adopted 
for $l in §2. With the single-exponential form, we set 
$ 2 (Ar,^) = exp(-|Ar|/T L ) for £ > L. An extension of 
this function to smaller scales gives, 



T(£) = T L . 

At smaller t, T(£) can be estimated using the veloc- 
ity scalings in the turbulent flow. For £ in the inertial 
range, we obtain T(£) by dividing £ by the amplitude of 

the turbulent velocity fluctuations at this scale, which is 
1/2 

(S a (£) + 2S nn {£)y / . Using the Kolmor gorov scaling, we 
have T{£) = C T e~ 1/3 ^ 2/3 , where C* T = (IICk/3)- 1 / 2 = 

— 1/2 

0.52C K with Ck being the Kolmorgorov constant. 
The factor, 11/3, is from the incompressibility relation 
S nn = 4S11/3 in the inertial range. Since Ck is ~ 2, we 
set Ct — 0.4. A similar value of Cr was adopted by 
Zaichik & Alipchenkov (2003). In the viscous range with 
£ <C 77, the flow velocity difference goes linearly £, and 
T{£) is expected to be constant. Lundgren (1981) pre- 
dicted that T{£) = V5t v , which was later confirmed by 
numerical simulations of Girimaji & Pope (1990). We 
thus take T(£) — \/5t v for £ -C 77 in our model. We 
will use the bridging formula for T(£) from Zaichik et al. 
(2006), 



$2 (At, £) — exp 



\Ar\_\ 
T{£))> 



(22) 



T(£) = T L 



where T{£) is essentially the correlation timescale or life- 
time of turbulent eddies of size £. For £ ^> L, we set 



l-exp(-(^0 32 (^) 



(£/ v y 



-1 -2/3 



(£/v7) 4 +(T L /(C T r,))6 



1/6 



(23) 
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which satisfies the scalings of T(£) in different scale 
ranges. 

One may also adopt a bi-exponential form for $2 (At, £) 
based on eq. ([5]) for the Lagrangian correlation function 
$ L (see §2). Replacing T L in eq. © by T{£) gives, 



* 2 (At,*) = 



1 



2yT - 2z 2 




V 7 ! - 2z 2 ) 



exp 



exp 



2|Ar| 



(1 + VI - 2z 2 )T(£) 



(1 - v 7 ! - 2z 2 



2|Ar| 



(1 - y/l - 2z 2 )T{£) 



(24) 



This bi-exponential form for $2(^1, £) was used in all 
the calculations in PP10. We will compute the predic- 
tions of the PP10 model using both the single- and bi- 
exponential correlation functions. We find the results 
from the two cases are close to each other, suggesting 
that the integral, eq. (|19p. is insensitive to the function 
form of $2 (At, £). After the integration, the dependence 
on $2 (At, I) is essentially condensed to a dependence 
on the correlation timescale T(£). This is similar to the 
case of the one-particle velocity, which is insensitive to 
the form of $i(Ar) (see §2). PP10 also considered the 
possible dependence of the parameter z on the length 
scale £. It was found that including a reasonable length 
scale dependence of z barely changes the model predic- 
tion. We will set z to be constant in this study. 

We next consider the simplified version of the PP10 
model represented by eqs. ([20| and (|2"Tj) . With the single- 
exponential form, eq. (|22l) . for $2, the response factor, 
cq. (f2i"j) , can be integrated, 



F(t) 



T(d) 



T{d) 



2r p T(d) 



(25) 

Since t is negative, the response factor, F(t), is domi- 
nated by the first term if T(d) 3> t p , and it approaches 
exp(T/T(d)) in that limit. On the other hand, for 

T(d) <C t p , the leading term is 2T ^ exp(T/T p ). Note 

that eq. ([231) does not diverge at T(d) — t p . Ap- 
plying the L'Hospital's rule shows that it converges to 
(| — ^)exp(T/T p ), as T(d) — > t p . Therefore, when 

numerically integrating eq. (|20|) . we set F(t) = (i — 
■^-) exp(r / t p ) for T(d) around t p . 

With the bi-exponential temporal correlation, eq. (|24|) . 
the response factor, F(t), can also be integrated analyt- 
ically. The integration is straightforward, but the result- 
ing function for F(t) is complicated and is thus omitted 
here. The predictions of the simplified model with single- 
and bi-exponential $2 (At, £) are also close to each other. 

3.2.2. Averaging over the direction of Ft 

We evaluate the angular average of Sij(R) over the 
direction of R. It follows from eq. (0 that Sij(R) = 
S aa (R)S ij + [Sn(R) - S nn (R)j RiRj/R 2 . The contraction 
of the tensor is Su(R) = Si\(R)+2S nn (R), and it does not 
have a direct dependence on the direction of R. There- 
fore, for the prediction of the 3D rms relative speed, we 



do not need to perform the angular average. However, 
for the radial and tangential components of the relative 
velocity, one must make an assumption for the direc- 
tion of R and compute the angular average of the term 

In PP10, we assumed that the direction of the sep- 
aration change, AR = R — r, caused by turbulent 
dispersion is completely random and isotropic. One 
can then insert R = AR + r into the above expres- 
sion of Sij(R) and take the average over the direction 
of AR. From the assumed isotropy of AR, we have 
(nARj) = and (ARiARj) ang = \{R 2 - r 2 )S lJ: and 
hence (R l R J ) ang ^ nr 3 + \(R 2 - r 2 )% (see PP10). The 
angular average (Sij(R)) ang is then given bj0, 



(Sij(R)) aIL g - $i; 



3i? 2 



Sn(R) 



3i? 2 



Snn(R) 



Sn{R) — S nn (R) 



R 2 ' 

(26) 



Clearly, the equation approaches SjAr) in the limit 
R — > r. PP10 showed that eq. (I26|) reproduces the 
Saffman- Turner formula for the radial and tangential 
relative speeds. In the limit with R 3> r, we have 
(£„-(#)) ang ~ § [Sn(R) + 2S nn (R)]. _ 

Here we make a simpler assumption than PP10: we 
take the direction of R (rather than AR) to be isotropic. 
This means (RiRj / R 2 ) ang = hSij, and we have, 



(27) 



oc 



(Sij(R)) aas = ~ [Sn(R) + 2S nn (R)} 

which suggests that the particle structure tensor S } 
Sij (see eq. (fl9^l). and hence (w 2 ) = (w 2 ) (see eq. ([TO])) for 
particles of any size. A comparison of the two assump- 
tions, eqs. ([26]) and (l27l) . shows that they differ only at 
R&r. 

As expected, the contraction, (Su{R)) ang , of both eq. 
flU and eq. (J27J), is equal to [S n (R) + 2S nn (R)], indicat- 
ing that the two assumptions give the same prediction 
for the 3D rms relative velocity. The only difference be- 
tween the two assumptions is the prediction for the radial 
and tangential components of the relative velocity in the 
small particle limit St <^L 1. In §3.2.4, we will compare 
the model predictions in this limit with both assump- 
tions. The angular average (Sij(d)) ang in the simplified 
model (eq. (|20p) can be evaluated in a similar way, and 
the resulting expressions are the same as eqs. (l26l) and 
(|27j) with d{r) replacing R(t,t'). 



3.2.3. The backward dispersion of particle pairs 

To complete the model, we finally specify the (rms) 
particle separation, d(r), as a function of t. The sep- 
aration of inertial particle pairs backward in time has 
not been explored in the literature. Fortunately, Bee 
et al. (2010) carried out a detailed numerical study of 
the forward-in-time pair dispersion of inertial particles. 

4 Rigorously, the amplitude, R, of R and hence S\\(R) and 
S n n{R) have a dependence on the direction of AR. However, the 
average of these quantities over the direction of AR is complicated 
and cannot be done analytically. For simplicity, we kept R, Sn and 
Snn fixed, and only accounted for the angular average of RiRj. 
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Following PP10, we use their results to guide the as- 
sumption for the backward dispersion. We first consider 
the separation behavior of inertial-range particles with 
Tj; ^ t p ^ Xl, and then discuss other particles at the end 
of this subsection. 

Bee et al. (2010) found that the separation of iner- 
tial particles shows different behaviors at early and late 
times. At early times, a clear ballistic phase is observed 
for particles with St <; 1. In this phase, the separa- 
tion increases linearly with time, and the phase lasts for 
about a friction timescale. The ballistic behavior is easy 
to understand: The particle velocity tends to be roughly 
constant for a memory timescale, r p . The reasoning also 
applies to the dispersion backward in time. We thus as- 
sume that, for particle pairs at an "initial" distance of 
r, the separation g?(t) in the time range — r p <^ r < is 
given by, 

d 2 { T )=r 2 + (w 2 )t 2 (28) 

where (w 2 ) is the 3D variance of the particle relative 
velocity at time 0. Note that the particle relative speed is 
actually what our model aims to predict. Therefore, the 
dependence of d(r) on (w 2 ) in the ballistic phase leads 
to an implicit equation for the particle relative velocity 
(see §3.2.4). 

Bee et al. (2010) also showed that, after a fric- 
tion timescale, the dispersion of inertial particles with 
T p Tl makes a transition to a tracer-like phase. In 
this phase, the particle distance is in the inertial range 
of the flow, and, like tracer particles, the separation vari- 
ance increases as time cubed, a behavior known as the 
Richardson's law. The Richardson behavior was found 
in the pair dispersion of tracer particles both forward 
and backward in time (Berg et al. 2006; see Appendix 
A) . It is thus likely to exist also in the backward separa- 
tion of inertial particles, and we expect it to connect to 
the ballistic phase at r ~ — r p . We therefore apply the 
Richardson's law for r <J — r p , 

d 2 (T)cge\r\ 3 (29) 

where g is called the Richardson constant and e is the 
average dissipation rate of the flow. Bee et al. (2010) 
did not report the value of g in the Richardson phase of 
inertial particle pair dispersion. As in PP10, we will take 
g as a parameter. We will measure g for the backward 
dispersion of tracer particles in our simulated flow (see 
Appendix A), and take the measured value as a reference. 
In our model, we will use a combined separation behavior 
that connects a ballistic phase and a Richardson phase 
at t — —Tp. 

The Richardson behavior would end when the separa- 
tion becomes much larger than the integral length scale, 
L, of the turbulent flow. At such a large distance, the 
flow velocities "seen" by the two particles is uncorrelated, 
and the separation of the particle pair is expected to be 
diffusive like in a random walk. It is thus appropriate to 
switch the Richardson behavior to a diffusive phase with 
d 2 (r) cx |t| at a separation d L. However, we find 
that the exact separation behavior at d L (or i? 3> L) 
does not affect the prediction of our model. This is be- 
cause at these large scales both the structure functions, 
S\\ and S\\, and the timescale, T{d) (or T(i?)), become 
independent of d (or R). Therefore, eq. (|T9| (or eq. (|20|) ) 
does not depend on the behavior of the separation once 



it becomes much larger than L. This is confirmed by 
the numerical solutions of eqs. (fl9|) and (|20|) . For conve- 
nience, we keep using the Richardson's law even after d 
exceeds L. 

The separation behavior discussed above is based on 
the simulation results of Bee et al. (2010) for particles in 
the inertial range. For simplicity, we will use the same 
behavior for all particles, although its validity is ques- 
tionable for small (r p <, t v ) and large (r p J> Tl) particles. 
For small particles with St <, 1, a ballistic phase is not 
clearly observed in the d 2 vs. time plots shown in Fig. 
5 of Bee et al. (2010). We expect that a short ballistic 
phase is likely to exist if one plots (d 2 — r 2 ) (instead of d 2 ) 
vs. time (see Fig. [I7]in Appendix A for the {d 2 — r 2 ) vs. 
time plot for tracer particle pairs). But the connection 
of this short ballistic phase to the Richardson behavior 
appears to be less clear than in the case of larger parti- 
cles (Fig. 5 of Bee et al. 2010). Fortunately, the relative 
velocity of St <C 1 particles is mainly contributed by the 
local flow velocity and does not have a strong dependence 
on the backward separation. With the assumed behav- 
ior, our model does give acceptable prediction for small 
particles. 

The problem of using the assumed behavior for par- 
ticles with r p ^> Tl is that the Richardson phase does 
not exit. The velocities of these large particles are un- 
correlated even at small distances (§3.1). Therefore, at 
timescales larger than r p , the separation is likely diffu- 
sive, i.e., d(r) cx |r|. Realistically, one needs to connect 
the ballistic phase to a diffusive behavior rather than a 
Richardson law at |r| ~ r p . However, it turns out that, 
at the end of the ballistic phase of these particles, the 
separation is already L. As discussed earlier, once 
the separation exceeds L, the exact separation behavior 
would not significantly affect the model prediction. This 
justifies using a combined separation behavior with a bal- 
listic and a Richardson phase also for r p 3> Tl particles. 

So far, in our assumption the initial distance, r, just 
provides a floor value for the separation d (eq. (|28|)). It 
is, however, possible that the value of r has additional 
effects on the separation behavior. Bee et al. (2010) only 
explored initial distances r above the Kolmorgorov scale, 
and it is not clear whether the separation behavior has a 
qualitative difference if r <> rj. In our model for the par- 
ticle collision speed, we are interested in the backward 
separation with r < 77, and it would thus be helpful to 
systematically investigate whether and how the separa- 
tion behavior changes as r decreases below r\. We defer 
such a study to a later work. In this paper, we assume 
that the two-phase behavior discussed above applies for 
any r. 

To constrain g in the Richardson phase, in Appendix 
A we measure g for the backward dispersion of tracer 
particle pairs in our simulated flow, which is used as a 
reference for inertial particles. The measured g for trac- 
ers in our flow at a limited resolution shows a dependence 
on r, suggesting that the Richardson constant for iner- 
tial particles may also depend on r. When comparing 
our model prediction to the simulation data for the par- 
ticle relative velocity at different r, we will adjust g to 
obtain best fits, and examine whether the best-fit values 
are consistent with the range of g measured from tracer 
particles. The Richardson constant for inertial particles 



11 



may also have a dependence on r p (or St) , which will be 
ignored for simplicity. 

3.2.4. Qualitative Behavior of Our Model Prediction 

Our model for the particle structure tensor, S v ij, is 
now complete. Here we discuss the qualitative behavior 
of our model prediction for the particle relative veloc- 
ity. We start by considering the 3D variance, (w 2 ). The 
contraction of eq. (flU)) gives, 

(w 2 )= — —[s n (R)+2S nn (R)]x 

J -oo T p J—oo T P 

$ 2 (r'~r,i?)exp (^j exp (j-j, (30) 

which is an implicit equation of (w 2 ) because R depends 
on (w 2 ) in the ballistic separation phase. In §6.1, we solve 
the equation numerically using an iteration method. 

The qualitative behavior of the model prediction for 
{w 2 } can be obtained by analyzing the integrand in eq. 
(|3"0|) . In the S-T limit (r p — > 0), the exponential cutoff 
terms, ^-exp(r/r p ) and exp(r'/r p ), in the integrand 

can be viewed as delta functions at r = and r 1 = 0, 
respectively. This suggests that (w 2 ) is approximately 
given by ~ (5 U + 2S* nn ) at R(0, 0). Since R(0, 0) = r, we 
have (w 2 ) = -^-r 2 for r in the dissipation rate, which is 
consistent with the 3D variance of the relative velocity 
from the S-T formula, eq. (jTTJ) . 

The analysis of eq. (|30j) for larger particles is more 
complicated. We first note that Su(R), S nn (R), and the 
timescale T(R) in the correlation function $2 are all in- 
creasing functions of R. Since R increases backward in 
time, the first factor in the integrand of eq. (|30|) increases 
with increasing |r| and |t'|. A larger T(R) also tends to 
increase the integral because $2 would allow contribu- 
tions from a broader range of time lag (At). Together 
with the exponential cutoffs, these suggest that the con- 
tribution to the integral peaks at r, r' ~ — r p . We de- 
note the particle separation at r = t' = — r p as R p 
(= R(— r p , — T p )), and refer to it as the primary distance. 
The primary distance provides useful insights to the pre- 
diction of our model. 

In the extreme limit of large particles with t p 3> Tl, 
R p is expected to be much larger that the integral scale, 
L, of the flow. We thus have S\\ = S nn — 2u' 2 , and 
T(Rp) = T L at Rp (see eqs. [3J and [23]). The expo- 
nential cutoff by $2 indicates that only the time pairs (t 
and t') that satisfy the constraint |t' — t| ~ Tl give a 
significant contribution to the integral. Since Tl <C t p , 
$2 reduces the range of r and t' that contributes to the 
integral by a factor of Tl/t p . Assuming the main con- 
tribution to the integral is from R ~ R p and accounting 
for the effect of <J> 2 , we find that (w 2 ) ~ 6w' 2 T l /t p . This 
is consistent with eq. (|12[) . meaning that our model cor- 
rectly reproduces the large particle limit. 

We next discuss inertial-range particles with 
T r] & T p <, Tl. For these particles, the primary 
distance, R p , corresponds to inertial-range scales of the 
turbulent flow. Using the Kolmorgorov scaling gives 

S\\,S nn oc i? 2 / 3 and T(R p ) oc R^/ 3 . From its defini- 
tion, Rp is roughly the particle distance at the time 
when the ballistic phase connects to the Richardson 



phase (see §3.2.3). We can thus assume that R p is 
determined by a ballistic separation of duration t p , i.e., 
Rp ~ (w 2 ) 1 / 2 T p . The effect of the $2 term depends on 
how T(Rp) compares to t p . If T(R p ) > t p , $ 2 is — 1 
for all time pairs in the range — t p <, t,t' ^ 0. On 
the other hand, if T(R P ) < t p , $ 2 provides a factor 
of T(R p )/t p , which follows from the argument used 
above for the large particle limit. We consider both 
cases and show that they actually lead to the same 
scaling of (w 2 ) with t p . In the first case, eq. ([30]) is 

approximated by (w 2 ) ~ [<Sn(i? p ) + S'nn(-Rp)] oc i? p ^ 3 . 

With R p ~ (w 2 ) 1 ^, we obtain (w 2 ) 1 / 2 oc t p 1/2 . 
In the second case with T(i? p ) < t p , we in- 
clude a factor of T(R p )/t p and estimate (w 2 ) as 

~ [Sn(Rp) + S nn (Rp)]T(Rp)/r p oc i? p /3 /r p . It is 
straightforward to see that setting R p ~ {w 2 ) x / 2 Tp in 

this estimate gives the same scaling, (w 2 ) 1 / 2 oc t p ^ 2 , as 
the first case. Therefore, whether T(R p ) is larger or 

smaller than t p , our model predicts a t p ^ 2 (or St 1 / 2 ) 
scaling for inertial-range particles. 

Using a similar argument, PP10 found that if the pri- 
mary distance is determined by the Richardson's law, 
Rp ~ (geTp ) 1 / 2 , the model also predicts a St 1 / 2 scaling 
for particles in the inertial range. Since both the bal- 
listic and Richardson behaviors yield a St 1 / 2 scaling, a 
combination of a ballistic and a Richardson phase pro- 
duces the same scaling (PP10). The St 1 ' 2 scaling has 
also been predicted by models of Volk et al. (1980), Cuzzi 
and Hogan (2003), Ormel and Cuzzi (2007), and Zaichik 
& Alipchenkov (2003). The derivation of the St 1 ' 2 scal- 
ing in all the models assumes a sufficiently broad inertial 
range in the turbulent flow. The scaling would not exist if 
the Reynolds number of the flow is low. In fact, the pre- 
dicted St 1 / 2 behavior has never been confirmed by simu- 
lations due to the low numerical resolution. PP10 showed 
that, to see the St 1 / 2 scaling, the Taylor Reynolds num- 
ber of the turbulent flow must be larger than ~ 300. This 
is higher than in the 512 3 simulation used in the present 
study, and thus the St 1 / 2 scaling is not observed. It 
appears likely that the existence of this scaling can be 
verified at a twice larger resolution. We will conduct a 
1024 3 simulation in a future work. 

The above analysis for the scaling behavior of (w 2 ) in 
different St ranges can be similarly applied to the sim- 
plified model, eqs. ([20]) and ([21"]) . The prediction of the 
simplified model is qualitatively the same as the original 
PP10 model. 

Finally, we examine the model prediction for the radial 
and tangential components of the relative velocity. The 
prediction for (w 2 ) and (w 2 ) depends on the angular av- 
erage of Stij in eq. (fl9]) (or eq. [20]). In §3.3.2, we made 
two assumptions, eqs. ([26]) and ([27]) . for the angular av- 
erage, and we discuss them separately. Inserting the first 
assumption, eq. (|26p , into eq. (|19p and comparing it with 
eq. ([T0| . we find, 
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(w 2 ) = 



dr 



dr' 



1 r 2 \ 
3 ^W 2 ) 



S U (R)- 



3 + 3i? 2 J nn 



(R) 



CXI) I — 



cxp — 



(31) 



In order to integrate the two equations, one needs to first 
solve eq. ([30j) for (w 2 ) due to the dependence of R on (w 2 ) 
in the ballistic phase. It is easy to show that, in the limit 



t p — > 0, i? — > r, and eq. (|3"Tj) reduces to 



S'il(r) and 

(w 2 ) = S nn (r), reproducing the Saffman- Turner formula, 



eq 



' 7) 1 



> r 

1 /„..2\ 



(ITT1) . For larger particles with r p ^> r, 
and thus eq. (|3"Tj) predicts that (i/; 2 ) = (w 2 / — g\uy /, 
Therefore, both (w 2 ) and (w 2 ) scale as S'i 1 / 2 for inertial- 
range particles and as St" 1 / 2 in the large particle limit. 

As discussed in §3.3.2, the second assumption, eq. ([27]) . 
for the angular average of Sry predicts that (w 2 ) = 
(w 2 ) = ■k(w 2 } for all particles. In the small particle limit, 



= J-r 2 , and thus 



eq. (f3T)f predicts ( 

^r 2 . This means that the prediction for the radial and 
tangential relative velocity of small particles by the sec- 
ond assumption differs from the S-T formula, although 
it reproduces the S-T prediction for the 3D rms. We will 
test the two assumptions against our simulation data. 

4. STATISTICS OF THE SIMULATED FLOW 

In this section, we describe the numerical method used 
in our simulation and discuss the statistical properties 
of the simulated flow. Our simulation was conducted 
in a periodic 512 3 box with a length of 2tt on each side. 
Using the Pencil code0 (Brandenburg & Dobler 2002, Jo- 
hansen, Andersen, & Brandenburg 2004), we evolved the 
hydrodynamic equations, 



dp 9 

m + dx~ {pUi) 



9uj 
dt 



9uj 
' dxj 



1 d 

pdxj 



pv 



9uj 



= 0, 

9uj 
9xi 



2 duk 
3 ij 9x k 



P OXt 



(32) 



with an isothermal equation of state, p = pC 2 . The 
sound speed is set to unity, i.e., C s = 1. The kinematic 
viscosity, v, is taken to be constant, v = 5 x 10 -5 . A 
large-scale force, fy, generated in Fourier space with wave 
numbers in the range 1 < k < 2 is applied to drive and 
maintain the turbulent flow. The driving length scale, 
Lf, is thus about 1/2 box size. The balance between 
the energy input by the driving force and the dissipation 
by viscosity leads to a statistical steady state with a ID 
rms velocity, u', of 0.05, or a 3D rms of 0.085. This 
weakly compressible flow is suitable for the application 
to turbulence in protoplanetary disks. At an rms Mach 
number of 0.085, the flow statistics is essentially the same 
as incompressible turbulence (Padoan et al. 2004, Pan & 
Scannapieco 2011). 

The integral length scale, L, in our simulated flow is 
found to be ~ 1, i.e., about 1/6 box size. It is about 
3 times smaller than the driving length scale, Lf. The 
integral scale, L, represents the (longitudinal) correlation 

5 http://pencil-code.nordita.org 



length of the velocity field, and we computed it from the 
energy spectrum, E(k), of the flow, using the relation 
L = 2^1 k' x E{k)dk (Monin & Yaglom 1975). The 
energy spectrum, E(k), is plot in the inset of Fig. Q. 
With L = 1, the large-eddy turnover time is T c ddy = 
L/u' = 20 in units in which the sound crossing time is 
2tt. 

The average energy dissipation rate per unit volume 
by the viscosity term is given by e = ^ (pv(diUj + djUi — 

^SijdkUk) 2 ), where p is the average density. In our 
weakly compressible flow, the density fluctuations and 
the velocity divergence can be neglected, and the aver- 



age dissipation rate can be estimated by e = v(lu ), where 
(oj 2 ) is the vorticity variance. We find that (w 2 ) = 0.92, 
implying that e ~ 4.6 x 10~ 5 . We also evaluated the 
dissipation rate from the 3rd-order longitudinal struc- 
ture function using Kolmorgorov's 4/5 law, (Au r (£) ) = 



for 



This latter method 

5 



1.04. The Kolmorgorov length scale is 
?? = (^ 3 /e) 1/4 = 0.0075, which corrc- 



, „ in the inertial range. 

gives a larger dissipation rate, e = 5 x 10~°, suggesting 
that a small fraction, ~ 8%, of kinetic energy is dissi- 
pated by numerical diffusion. The effective viscosity is 
thus larger than the adopted value by the same amount. 
We take the effective viscosity to be 5.4 x 10~ 5 and use it 
in our estimates of the Kolmorgorov scales. We compute 
the Kolmorgorov timescale from the vorticity variance as 

T V = (W 2 )- 1 / 2 = 

estimated to be 

sponds to ~ 0.6 cell size of the computation grid. The 
Kolmorgorov velocity scale is u v = {ye) 1 /* = 0.0072 in 
units of the sound speed. 

The Reynolds number of our simulated flow is Re = 
u'L/v ~ 1000. A more commonly-used Reynolds num- 
ber in turbulence studies is the Taylor Reynolds number, 
Re\ = u'X/v, based on the Taylor micro length scale, A. 
The Taylor length scale is defined as A = (15u' 2 / (w 2 )) 1 / 2 . 
We find that A = 0.2 in our simulated flow, and thus the 
Taylor Reynolds number Re\ ~ 200. From their defini- 
tions, we have v! /u v = (Rex/VTE) 1 / 2 , which is useful for 
normalizations . 

4.1. The Lagrangian Correlation Function and the 
Timescales 

To study the Lagrangian statistics, we integrated the 
trajectories of 33.6 million tracer particles with zero in- 
ertia in the simulated flow. The total number of tracer 
particles corresponds to an average number density of 1 
particle per 4 computational cells. To obtain the particle 
velocity inside a cell, we selected the triangular-shaped- 
cloud interpolation method already implemented in the 
Pencil code (Johansen and Youdin 2007). We output 
the particle positions to a data file in each time interval 
of O.lr,,. The Lagrangian correlation function, $l(At), 
is computed as the average of the velocity correlation, 
(u I (X(t),t)M l (X(t + At),* + Ar))/3u' 2 , along the tra- 
jectories, X(t), of all particles. We considered both pos- 
itive and negative time lags At, corresponding to La- 
grangian trajectories forward and backward in time, re- 
spectively. Our data confirmed that is an even func- 
tion of At, as expected from statistical stationarity (see 
§2). We find that the Lagrangian correlation timescale, 
T L (= / $ l (At)(2At), is ~ 15, which is about 0.75 eddy 
turnover time, T e ddy This is consistent with the Simula- 
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Fig. 2. — Lagrangian (<1?l; circles) and Eulerian ("^e; triangles) 
temporal correlation functions in our simulated flow. The time lag, 
At, is normalized to the Kolmorgorov timescale, t v . The solid line 
shows the best fit for <E>l using the bi-exponential function, eq. 

tion result of Yeung et al. (2006). Since t v = 1.04 in our 
simulated flow, we have Tl = 14. At v . 

The Lagrangian correlation, $l, in our flow is plot as 
circles in Fig. [21 where the time lag, |Ar|, is normalized 
to the Kolmorgorov timescale, t„. The solid line shows 
the bi-exponential function, eq. ([5]), given in §2. The pa- 
rameter z is set to 0.3, which suggests that the Taylor 
micro timescale, tt, is ~ 4.3r, ; . This value of tt cor- 
responds to an acceleration variance, a 2 ~ 5.2(u i? /t 7; ) 2 . 
The bi-exponential function matches very well the sim- 
ulation data. On the other hand, we find that a single 
exponential function could not give a satisfactory fit to 
■i'l, 

We also considered the Eulerian temporal correlation 
function, $e(At), in the simulated flow. It is computed 
as the average, (ui(x, t)ui(x, t + At)) /3m' 2 , over all grid 
points a;. The result is plot as triangles in Fig. [2j $e(At) 
is also symmetric about At = 0. As seen in Fig. [2J S'e is 
smaller than the Lagrangian correlation $l at small time 
lags, and then becomes larger at |Ar| ^ 8r I( ~ 0.5TL. 
Due to the slower decrease of <!>e at large time lags, 
the Eulerian correlation time, Te = J §e(At)cIAt, is 
slightly (10%) larger than T L . We find that T E = 15.%. 
The Eulerian correlation function is of interest for large 
inertial particles with t p ^> Tl. Due to their large in- 
ertia, these particles have small velocities and thus may 
stay around as the flow sweeps by. Therefore, unlike 
small particles, the temporal series of the flow velocity 
"seen" by the large particles may be better described by 
the Eulerian velocity. This suggests that, for r p > Tl, it 
may be appropriate to replace the Lagrangian correlation 
used in our model by the Eulerian correlation. However, 
the Eulerian correlation function and timescale are quite 
close to the Lagrangian ones, and using the Lagrangian 
correlation for all particles in our model gives satisfactory 
predictions for the 1-particle velocity and the 2-particle 
relative velocity at any value of St (§5 and §6.1). 

We give a brief summary of the relevant timescales 
in the simulated flow. We list them in an increasing or- 
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Fig. 3. — Longitudinal (open circles) and transverse (filled circles) 
structure functions measured in our simulated flow. The structure 
functions are compensated by the Kolmorgorov scaling in the iner- 
tial r ang e. The solid and dashed lines are fitting functions, eqs. 
and i'S'St . for S\\ and S n n, respectively. The inset shows the energy 
spectrum of the simulated flow, compensated by e 2 / 3 fc -5 / 3 . 

der. The smallest timescale is the Kolmorgorov timescale 
r^, and we use it as a reference timescale. The Tay- 
lor micro scale, Tt, was found to be 4.3^ from the bi- 
exponential fit for the Lagrangian correlation function. 
The next timescale is the Lagrangian correlation time, 
Tl, which is 14. 4t v . The Eulerian correlation time is 
slightly larger, Te ~ 15.9t, ; . The large eddy turnover 
time, T ec idy, was measured to be ~ 19.2^. Another 
commonly- used timescale is the dynamical time, Td yn , 
defined as the forcing length scale, L[ , divided by the 3D 
rms velocity (\/3u'). We find that Td yn = 35^. 

4.2. The Flow Structure Functions and Energy 
Spectrum 

In Fig. [31 we show the longitudinal (Sn; open circles) 
and transverse (S nn ; filled circles) structure functions in 
our simulated flow. The structure functions are mea- 
sured from the velocity differences along the 3 direc- 
tions, e\, e2 and e 3 , of the simulation grid. For S\\(£), 
we computed and averaged the variances of Aun(£)(= 
ui(x + lei) — ui(x)), Au22{£) and Ait33(£) over all the 
grid points, x. Similarly, S nll (£) is obtained by averaging 
the variances of Aui2(£)(= ui(x+£e2) — ui(x)), Aiti3(£), 
Au2i(£), Au 23 (£), Au 31 (£) and Au 32 (^ 

As discussed earlier, Kolmorgorov's similarity theory 
predicts that S\\(£) ~ Ck(e^) 2 / 3 for £ in the inertial 
range. We thus compensated the structure functions by 
(e£) 2 / 3 in Fig. [3] A limited inertial range is seen in both 
Sn and S nn . Clearly, the Kolmorgorov constant Ck is 
about 2. In this range, the scaling exponent for S\i is 
found to be slightly larger than 2/3, while the slope of 
S'nn is close to 2/3. The ratio of the two structure func- 
tions in the inertial range is about 1.25, slightly smaller 
than the value, 4/3, expected from the incompressibil- 
ity condition (see §3). This is perhaps because our flow 
is weakly compressible. Another possible reason is that 
the inertial range is too short to allow an accurate mea- 
surement of this ratio. Both structure functions become 
smooth, i.e., cx £ 2 , as £ decreases toward the Kolmorgorov 
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scale, and approach 2m' 2 in the limit I >• L (§3). 

The solid line in Fig. [3] is the connecting formula, eq. 
©, for Sn (§3). We set C K to 2 in the formula. The line 
gives a fairly good fit to the data points. As discussed in 
§3, with the connecting formula for Sn, one may obtain a 
fitting function for S nn using the incompressibility rela- 
tion S nn = S\\ + T;£dS\\/d£. However, the fitting function 
obtained this way overestimates the data points for S nn 
in the inertial range, perhaps because the incompress- 
ibility condition does not exactly hold in our flow (see 
above) . For a more accurate fitting function, we adopted 
a separate connecting formula for S nn , 
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(33) 



where Ckn is the scaling coefficient for S nn in the inertial 
range. This connecting formula correctly reproduces the 
scaling behaviors of S nn in different scale ranges. Its form 
is slightly different from eq. (JU) for Sn. The dotted line in 
Fig. [3] corresponds to eq. (|33|) with Ckn — 1.25Ck = 2.5. 
We will use eqs. © and (|3"3")l in the computation of our 
model prediction for the particle relative velocity. 

The inset of Fig. [3] show the energy spectrum, E(k), 
of our flow. The Kolmorgorov theory predicts E{k) = 
i"Ce 2 / 3 fc~ 5 / 3 in the inertial range, and we compensated 
the spectrum by e 2 / 3 fc -5 / 3 . The power-law range (3 < 
k < 10) in the spectrum appears to be shorter than in 
the structure functions. The constant K is measured to 
be ~ 1.7, consistent with previous studies (Ishihara et al. 
2009). It is also consistent with the relation K = 0.76Ck 
(Monin & Yaglon 1975), as the Kolmorgorov constant, 
Ck, for Sn was found to be ~ 2. 

5. ONE-PARTICLE ROOT-MEAN-SQUARE VELOCITY 

In our simulation, we included 14 species of inertial 
particles of different sizes. The friction timescale of the 
particles spans about four decades from ~ O.lr^ to 54TL 
(~ 800r^), covering the entire scale range of the simu- 
lated flow. The friction timescale is equally spaced, in- 
creasing by a factor of two in each successive species. The 
number of particles contained in each species is 33.6 mil- 
lion, which is the same as the tracer particles used in the 
study of Lagrangian statistics (§4.1), corresponding to an 
average particle density of one per 4 computational cells. 
Similar to the integration of the tracer particle trajecto- 
ries, when evolving the equation of motion (eq. [T]) for in- 
ertial particles, we adopted the triangular-shaped-cloud 
method to interpolate the flow velocity. Our simulation 
run lasted for about 35Tl, which is much larger than 
the friction timescale of all particles except the largest 
ones with r p ~ 27 and 54TL. The dynamics of all par- 
ticles in the first 12 species (r p <, 14Tl) is thus fully 
relaxed by the end of the simulation, and we measured 
the statistics of the particle velocity using the last few 
snapshots covering 5 — 6 Tl. On the other hand, it is 
uncertain whether the largest particles from the last two 
species are well relaxed. From the analysis of our sim- 
ulation data, we find that the velocity statistics of the 
largest particles do reach a steady state toward the end 
of the run, suggesting that these particles may also be 




Fig. 4. — One-particle rms velocity, t)', as a function of the Stokes 
number, St. The data points show the simulation data. The dotted 
and solid lines are the predictions, eq. $Q and eq. 10 , of our model 
using single- and bi- exponential forms for the temporal correlation 
function, respectively. In both cases, Tl is set to 15 At v , and the 
parameter z in the bi-exponential case is set to 0.3. 

dynamically relaxed. 

In Fig. 21 we show our simulation result for the ID rms 
particle velocity, v' , as a function of the Stokes number, 
St. We normalized v' by the flow rms velocity, u' . The 
dotted and solid lines are predictions, eq. ((4} and eq. 
© , of our model using single- and bi-exponential forms 
for the temporal correlation function, respectively. As 
discussed in §2, the model approximates the trajectory 
correlation function, $i, by the Lagrangian correlation 
function, $l- For the bi-exponential case, we set the pa- 
rameter z — 0.3, which best fits <!>l measured from the 
Lagrangian trajectories of tracer particles (see §4.1). The 
two lines almost coincide, indicating that the model pre- 
diction for v' is insensitive to the exact form of the corre- 
lation function, and depends essentially only on the cor- 
relation timescaleQ- In both curves, we set Tl = 15At v , 
generally consistent with the directly measured value of 
14.4t, ; (see §4.1). 

As discussed in §4.1, the flow velocity "seen" by the 
large particles with t p J> Tl may be closer to Eulerian 
than Lagrangian, and thus using <I>l for the trajectory 
correlation <I>i is physically not well justified (§2). How- 
ever, the assumption is validated by our simulation data 
for all particles. This is because, first, based on our 
model prediction, v' is controlled mainly by the corre- 
lation time, not by the exact shape of the correlation 
function, and, second, the Eulerian correlation timescale, 
Te, is close to Tl. Therefore, whether $i is better de- 
scribed by the Eulerian correlation function $e or <E>l, 
the 1-particle rms velocity would be essentially the same, 
justifying using <1>l for all particles. We point out that 

7 Unlike the case of the 1-particle velocity, the choice of 3>i is 
crucial for predicting the relative velocity between different par- 
ticles in the S-T limit (PP10). In that case, adopting the bi- 
exponential form of <I>l for $i is needed to reproduce the accelera- 
tion contribution to the relative velocity, while a single exponential 
form cannot correctly capture the effect of the flow acceleration on 
4>L at small time lags. 
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the best-fit correlation timescale, 15. At v , in Fig. 0] is ac- 
tually in between the measured values of TL (14.4r^) and 
Te (15.9T,,). This suggests that the temporal correlation 
of the flow velocity "seen" by a large particle is likely 
in between that along a Lagrangian trajectory and that 
at a fixed Eulerian point. It is interesting to note that, 
even though a single-exponential function does not fit 
well either $e or $l (§4.1), our model prediction with 
an exponential correlation is in good agreement with the 
simulation result for the 1-particle rms velocity. 

We also computed the probability distribution function 
of the 1-particle velocity. The velocity of small particles 
is expected to be Gaussian because it simply samples 
the 1-point PDF of the flow velocity, which is close to 
Gaussian. The velocity of large particles with r p ^> Tl 
would also be Gaussian because their equation of motion 
is essentially a Langevin equation. We find the 1-particle 
velocity PDF is actually nearly Gaussian at all St. The 
far PDF tails are slightly thinner than Gaussian probably 
due to the sub-Gaussian tails of the flow velocity and/or 
the insufficient sample size at high tails. 

6. THE RELATIVE VELOCITY OF INERTIAL PARTICLES 

We finally explore the statistics of the relative veloc- 
ity of inertial particles in our simulation, focusing on 
the monodisperse case with identical particles. For each 
species (St), we measure the relative velocity of particle 
pairs at three given distances, r = lr), 0.5?7 and O.2577. A 
study of smaller values of r is desirable, as the size of dust 
particles in protoplanetary disks is much smaller than the 
Kolmorgorov scale. However, at smaller r, the number of 
particle pairs available in our simulation becomes more 
limited and does not allow accurate statistical analysis. 
For each St and r, we search the simulation box for all 
particle pairs in a distance shell from r—5r to r+Sr, with 
Sr set to O.Olr. We decompose the relative velocity, w, 
of each selected particle pair into radial and tangential 
components, and then analyze their statistics separately. 

To compute the radial and tangential components, we 
set up a local coordinate system (e[, e' 2 and e' 3 ) for 
each pair. The direction e\ is chosen to coincide with 
the separation, r, of the two particles. In terms of the 
unit base vectors, e\, e 2 and e 3 , of the simulation grid, 
e\ is expressed as cos#cos</>ei + cos # sin <^e2 + singes, 
where sin6> = r^/r, cos 9 = (r 2 + r\) x l 2 jr, cos0 = 
T\j(r\ + r|) 1,/2 , and sin0 = r 2 j(r\ + r^) 1 / 2 . The ra- 
dial relative velocity is calculated as w T — w ■ e[ . For the 
two tangential directions, we set e' 2 = — sin §e.\ +cos <f>e 2 
and e 3 = — sin 6* cos <j>e\ — sin#sin</>e2 + cos#e3, which 
are obtained by two consecutive rotations of the orig- 
inal coordinates. The first rotation is about e 3 by <f>, 
which moves e 2 to e' 2 , and the second rotation is about 
e' 2 by —9, which further brings the original base vectors 
ei and e 3 to e[ and e 3 . We calculate the tangential 
components of the relative velocity w as Wt2 = w ■ e' 2 
and u>t3 = w ■ e 3 . The PDFs of u>t 2 and Wts are found 
to be almost the same, as expected from the statistical 
isotropjQ- We thus take the PDF of a tangential com- 
ponent, wt, to be the average PDF of Wt2 and Wt3- The 

9 When selecting the local coordinate system, one may also per- 
form a third rotation about e'j by an arbitrary angle. This changes 
e' 2 and e^. However, from the statistical isotropy, Wt2 and Wt3 
would be statistically invariant with respect to this third rotation. 



variance of w t is calculated as (w 2 ) = | ((w 2 2 ) + (w 2 3 )) . 
We discuss the measured statistics of w T , w t and the 3D 
amplitude |io| in the following sections. 

6.1. The Root-mean-square Relative Speed 
In Fig. [5j we show the simulation result for the 3D 
rms, (w 2 ) 1 / 2 , of the relative velocity as a function of 
the Stokes number. The data points correspond to the 
measured relative velocity for particle pairs at a distance 
of li]. We normalized the relative speed by the Kol- 
morgorov velocity, u<q. One may easily convert the nor- 
malization to the rms flow velocity u' using the relation 
\j 2 

v! = (Re\/V~L5) u v = 7ur). At small St, the 3D rms 
relative speed is roughly constant, and its value is con- 
sistent with the prediction, (w 2 ) 1 / 2 = Ur,/V3, of the S-T 
formula. The relative speed starts to rise at St ~ 1, as 
the contribution from the backward separation of parti- 
cle pairs becomes important. For the largest particles, 
we find that (w 2 ) 1 ' 2 ~ vW '(TL/Tp) 1 / 2 , with T L ~ 14r^, 
in agreement with eq. (TH21) for the large particle limit, 
r p 3> TL . Like the earlier result for the 1-particle rms ve- 
locity, this provides a validation for using the Lagrangian 
correlation time for the trajectory correlation, $1, of 
large particles with r p 3> Tl, even though their trajec- 
tories may significantly deviate from Lagrangian tracers. 
The predicted St 1 / 2 scaling for inertial-range particles by 
various models is not observed due to the limited inertial 
range of the simulated flow. 

The solid curve in Fig. §5§ is the prediction of the PP10 
model, and it is obtained by numerically solving eq. (|30p 
with an iteration method. We computed the trajectory 
structure tensor, Sry, using eq. ((8|) for Si\, eq. (|33|) for 
S^n, and eq. (|23]l for T(£), respectively. The parameters 
used in these equations are Ck = 2, Cxn = 2.5, Ct = 0.4 
and Tl = 14.4t, ; . A bi-exponential form, eq. ([M]). is 
adopted for the temporal correlation <J> 2 - We used a 
two-phase behavior for the particle separation backward 
in time (§3.4.3), connecting the ballistic and Richard- 
son phases at r c = — r p . We set d 2 (r) — r 2 + (w 2 )t 2 
for t > t c , and then switch to the Richardson's law, 
d 2 (r) = d 2 (r c ) + ge(\r\ 3 - |t c | 3 ), for r < t c , where 
d 2 (r c ) — r 2 + (w 2 )t 2 . To fit the simulation data, the 
Richardson constant, g, is set to 1.6. The solid line is 
in good agreement with the data points, confirming the 
validity of the physical picture of our model. Adopting a 
larger value for g could further improve the fitting quality 
at intermediate St. 

In Appendix A, we investigated the backward separa- 
tion of tracer particles in our simulated flow, and find 
that g is in the range 0.5 <, g < 1.2. Therefore, the g 
value used in the solid line in Fig. [5] is significantly larger 
than that for tracer particles. It is possible that the 
backward separation of inertial particles in the Richard- 
son phase is indeed faster than tracers, which can be 
tested by a direct numerical study of inertial particle 
dispersion. Another possibility is that the accuracy of 
our model for the trajectory structure tensor, Stvj, is in- 
adequate. First, there are uncertainties in our estimate 
for the correlation timescale, T(£), in $ 2 . We associated 
the timescale with the eddy turnover time at the typical 
particle separation, R(t,t'), between r and r'. This ap- 
proximation is likely not of high quantitative accuracy. 
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Fig. 5. — The 3D rms relative velocity, (ui 2 ) 1 / 2 , as a function of 
the Stokes number, St. The data points show the simul atio n result. 
The solid line is the prediction of the PPIO model (eg. 1301) using a 
bi-exponential temporal correlation $2 and a two-phase separation 
behavior with g = 1.6 in the Richardson phase. The dotted line 
corresponds to the same model but with a single-exponential <3?2- 
The line is barely visible because it almost coincides with the solid 
line. The dot-dashed line assumes ballistic separation at all ti mes . 
The dash ed line is the prediction of the simplified model, eqs. I|20ll 
and 1 1211) . using single-exponential $2 and a two-phase separation 
behavior with g = 1.0. 

At inertial-range scales, we set T(t) = C^l'^^i 2 ^ with 
Ct — 0.4. There may be an order-of-unity uncertainty in 
the estimated value of Ct- Our model prediction for in- 
termediate St increases with Ct- If a larger value of Ct 
were adopted, we would obtain a good fit to the simula- 
tion data with a smaller value of g. Another uncertainty 
in our model for Stij is that we neglected the correlation 
between the particle distance, R, and the fluctuations in 
the flow velocity difference, Au(R), seen by the particle 
pair (see §3.2). This assumption tends to underestimate 
the relative velocity. Accounting for the correlation may 
lead to a better agreement of the model prediction with 
the data using a smaller value of g. In §6.3, we will show 
evidence for this correlation from our simulation result 
for the PDF of the particle relative velocity. 

The dotted line in Fig. ((5) is the prediction of the 
same PPIO model, but with a single-exponential form 
(eq. (|22])) for $2- The dotted line is actually not dis- 
tinguishable from the solid line, confirming the earlier 
statement in §3.2.1 that our model prediction for the rel- 
ative speed is insensitive to the function form of $2- This 
leaves us some freedom for the choice of $2, as long as 
the correlation timescale is properly evaluated. In par- 
ticular, it provides a justification for approximating $2 
by the same function form, e.g., the form of $l, for all 
particles, although realistically the form of $2 may have 
a dependence on the particle inertia. 

The dot-dashed curve plots the prediction of the PPIO 
model assuming that the particle separation is ballistic 
with d 2 (r) = r 2 + (w 2 )t 2 at all times. The model is 
otherwise the same as the solid line. A pure ballistic 
separation is not realistic, and we show it here just to 
illustrate whether the Richardson separation phase pro- 
vides important contribution to the relative velocity. At 



0.5 St <z 5, the dot-dashed line significantly under- 
estimates (u> 2 ), and, from St ~ 5, it becomes close to 
both the data points and the solid line using a two-phase 
separation behavior. A possible explanation for this is 
that, for 0.5 St 5S 5, the relative velocity receives a 
significant contribution from the Richardson phase, even 
though this phase occurs at times beyond the particle 
memory timescale, i.e., at r < — r p . In that case, ac- 
counting for this phase would then be necessary for par- 
ticles with intermediate St. However, the validity of the 
above interpretation is subject to future tests. There 
is the possibility that the discrepancy between the dot- 
dashed line and the data points is mainly caused by var- 
ious uncertainties in our model for the trajectory struc- 
ture tensor, S^ij (see above). It is possible that, with 
an improved model for Stij, a pure ballistic separation 
may give a prediction that better matches the simulation 
data for St in the range 0.5 & St ^ 5. 

The dashed line is the prediction of the simplified 
model, eqs. (f2T))) and (|2"B"j) . using a single-exponential $2- 
As before, the model with a bi-exponential $2 gives al- 
most the same prediction. The same two-phase separation 
behavior as in the solid line for the original PP10 model 
is adopted. For the simplified model, the best-fit value of 
g is found to be ~ 1, which is close to the g values mea- 
sured from tracer particles. The simplified model also 
fits the data better for intermediate St, although the as- 
sumption made in the model is physically not better than 
the original PP10 model. The simplified model may be 
a preferred choice, as its prediction is easier to compute. 

We find that the Stokes number, St m , at which the 
rms relative velocity peaks is ~ 30, corresponding to a 
friction timescale twice larger than Tl. The peak value 
of the 3D rms relative velocity is ~ 6.2^, twice smaller 
than the 3D rms velocity (\/3u') of the flow. We give 
an explanation for the behavior of the peak relative ve- 
locity using the qualitative analysis of our model pre- 
diction given in §3.2.4. The analysis was based on the 
primary distance R p , estimated as R p = (w 2 ) 1 / 2 r p . R p 
generally increases with t p . Around the peak of the rel- 
ative velocity, r p ~ 30r„ and (w 2 ^ 1 ^ 2 



R p ~ 2OO77. From eq. 



6.2U,;, and thus 



the correlation time, T(£), 



at £ ~ 2OO77 is about 14^, which is close to Tl. For 
St ;> 30, T(Rp) would be constant and ~ Tl. Con- 
sequently, the $2 term in eq. (f30"|) provides a factor of 
Tl/ t p for all particles with St ^ 30 (see §3.2.4). Using 
the same analysis as in §3.2.4, one can show that this 
factor causes the relative speed to decrease with t p , even 
though the structure functions S\i(R p ) and S nn (Rp) are 
still increasing with R p for R p ~ 200ry (see Fig. [3| . For 
particles with St ^ 30, both the structure functions and 
T(Rp) decrease with decreasing R p , and thus the rela- 
tive speed would decrease with decreasing r p . Therefore 
a peak forms at St m ~ 30. We find that, for particles 
with St ~ 30, the amplitude of the flow velocity differ- 
ence at the primary distance (i? p ~ 2OO77) is smaller than 

^/3u' , and this is responsible for why the maximum rel- 
ative velocity is significantly lower than the rms velocity 
of the flow. The discussion here shows that the relative 
velocity is the largest for the particles whose primary dis- 
tance R p corresponds to the size of turbulent eddies with 
lifetime ~ Tl. Clearly, the backward particle separation 
plays an important role in determining the peak Stokes 
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number, St m . 

We give a brief discussion on the model of Volk et 
al. (1980) and its later developments and show that it is 
likely that their prediction does not well match our simu- 
lation data. In these models, the predicted relative speed 
reaches its maximum value when r p is equal to a large 
eddy time, tj, (e.g., Markiewicz, Mizuno & Volk 1991, 
Cuzzi and Hogan 2003, Ormel & Cuzzi 2007). The defini- 
tion of £l in these studies is different from the timescales 
used in the current work, and it is not clear whether, 
using parameters appropriate for our simulated flow, the 
model may correctly produce a peak at St m = 30. One 
may need to tune the value of t^ in the model such that 
the peak of the relative speed occurs at the right Stokes 
number. Another issue is that, in both the original Volk 
et al. model and its later extensions, the peak relative 
speed is predicted to be equal or close to the rms veloc- 
ity of the flow. This means that the model overestimates 
the relative speed by a factor of 2 around the peak, as- 
suming the peak position, St m , is correctly produced. A 
physical problem of the Volk et al. model has been dis- 
cussed in the Introduction, and the interested reader is 
referred to PP10 for details. One possibility is that the 
performance of the Volk et al. (1980) model may improve 
as the Reynolds number of the flow increases. This will 
be tested in a future work with higher-resolution simula- 
tions. 

6.1.1. Dependence on the Particle Distance 

Fig. [5] plots the 3D rms relative velocity for particle 
pairs at different distances, r. The squares, circles and 
diamonds correspond to r = 1, 0.5 and 0.25?7, respec- 
tively. The relative velocity shows shows a r-dependence 
at St < 6, while it is independent of r for St > 6 particles. 
In the context of our physical picture, this is because the 
friction time, r p , of St ^> 6 particles is long enough that 
the backward particle separation after a duration of ~ r p 
is insensitive to the "initial" value, r. On the other hand, 
the relative speed of smaller particles relies on the flow 
velocity difference they saw in the near past, when the 
particle separation was still dependent on r. 

The solid and dotted lines are predictions of the PP10 
model with bi-exponential correlation function $2 and 
the simplified version with single-exponential $2, respec- 
tively. The lines for r = lrj have already been shown in 
Fig. ([3]), and the Richardson constant, g, was set to 1.6 
and 1.0, respectively, in the two models. We find that, 
at smaller r, the best-fit value of g becomes smaller. For 
the PP10 model, we adopted g = 1.3 and 1.0 for r = 0.5 
and 0.2577, respectively. The decrease of g with decreas- 
ing r is consistent with our result in Appendix A for the 
pair separation of tracer particles. The backward disper- 
sion of tracer pairs was found to be slower for smaller 
r. Similarly, the g value used in the simplified model 
also decreases with decreasing r. In the dashed lines for 
r = 1, 0.5 and 0.25?7, the value of g is set to 1, 0.7 and 
0.5, respectively. 

The relative speed of the smallest particles (St ~ 0.1) 
in our simulation appears to be larger than the second 
smallest ones (St ~ 0.2), especially at smaller r. Slight 
dips are thus seen at St ~ 0.2 (Fig.[BJ). This is in contrast 
to the S-T formula, which predicts the relative speed 
at a given r is constant at sufficiently small St. These 
dips are not expected from the physical picture of our 
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Fig. 6. — The 3D rms relative velocity, (w 2 ) 1 / 2 , at at r = 1 
(squares), 0.5 (circles) and 0.25?7 (diamonds). Solid and dashed 
lines are predictions of the PP10 model with bi-exponential $2 and 
the simplified version with single-exponential $2, respectively. The 
two-phase separation behavior connecting at r = — t p is adopted 
in both models. In the PP10 model, the Richardson constant g 
is set to 1.6, 1.3 and 1.0 for solid lines from top to bottom. The 
corresponding g values used in the simplified model are 1, 0.7 and 
0.5, respectively. 

model either, and their existence is thus questionable. 
One possibility is that the rise of the relative speed to- 
ward St ~ 0.1 is a numerical artifact. This suspect is 
based on the consideration that the trajectory integra- 
tion in our simulation is likely less accurate for smaller 
particles. The accuracy of the trajectory computation 
depends on the integration time step relative to the fric- 
tion timescale, t p . The time step is the same for all 
particles in our simulation, so the accuracy is expected 
to be lower for smaller particles. Since we are not aware 
of a convincing physical mechanism that may cause a dip 
at St ~ 0.2, we tentatively assume the rise of the relative 
speed at St ~ 0.1 is due to a numerical artifact. This 
hypothesis will have to be tested with future simulations 
with a better temporal resolution of the trajectories of 
small particles. 

As mentioned earlier, the 3D rms relative velocity mea- 
sured in our simulation for small particles at r = I77 is 
consistent with the S-T prediction. The S-T formula also 
predicts that (w 2 ) 1 / 2 scales linearly with r in the St -C 1 
limit. This linear scaling is, however, not confirmed by 
the data for the smallest particles in the simulation. A 
rough power fit for the rms relative speed as a function 
of r gives (w 2 ) 1 / 2 cx r 78 at St = 0.1 - 0.2. This means 
that, for particles with St ~ 0.1 — 0.2, the S-T formula 
is already invalid at r <, 0.577. At a given St, a criti- 
cal particle distance is expected, below which the linear 
scaling does not apply. The physical reason is that, as r 
decreases, the local flow velocity difference across r be- 
comes smaller, and it is easier for the particle memory of 
the flow velocity difference in the past to provide a sig- 
nificant contribution, which tends to invalidate the S-T 
prediction. Equivalently, at a given r, the S-T formula 
is valid only below a critical St. In Fig. [6l the lines for 
r <, O.577 show that, as St decreases to ~ 0.2, the relative 
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speed predicted by our model is not flat yet, suggesting 
a significant contribution from the particle memory. For 
r <, O.577, the S-T prediction may be recovered at St 
much smaller than 0.2, assuming that the rise of the rel- 
ative speed toward St ~ 0.1 in the data points is indeed 
a numerical artifact. 

Our model does not directly account for the sling ef- 
fect (Falkovich et al. 2002, Falkovich and Pumir 2007) 
or the related caustic formation (Wilkinson & Mehlig 
2005; Wilkinson et al. 2006; Gustavsson k Mehlig 2011). 
The effect of slings is usually considered in the context 
of small particles with St <, 1. As mentioned in the In- 
troduction, the effect corresponds to crossing of particle 
trajectory that occurs at fluid streamlines with high cur- 
vature or local flow regions with large velocity gradient. 
The trajectory crossing contributes to the relative veloc- 
ity of particles especially at small distances. In our physi- 
cal picture, the effect of slings or caustics could be viewed 
as a contribution to the backward particle separation, 
and, in principle, it may be incorporated in our model 
by quantifying the frequency and contribution of such 
events. Based on the model of Wilkinson et al. (2006) 
and Gustavsson k Mehlig (2011), at a given small St, the 
effects of caustics would dominate the S-T contribution 
at a sufficiently small distance, r. Falkovich and Pumir 
(2007) showed that the sling effect is already significant 
at St ~ 0.2. In Fig. [6j we see that our model prediction 
underestimates the relative velocity of St = 0.2 particles 
at r <, O.S/;. This is a sign of the contribution from slings 
or caustics. We will discuss the effect of caustics on the 
particle collision rate in §6.2. 

6.1.2. The Radial and Tangential Relative Speeds 

In Fig. [71 we show the rms relative speeds in the radial 
(circles) and tangential (diamonds) directions for particle 
pairs at r = 1, 0.5 and 0.257/. At St <, 1, the tangential 
rms speed, (u; 2 ) 1 / 2 , is slightly larger (by ~ 10%) than 
the radial rms, (w 2 ) 1 ^ 2 . This difference is significantly 
smaller than the prediction of the S-T formula, eq. (|TT|) . 
that predicts that in the St <C 1 limit the tangential rms 
relative speed should be larger than the radial rms by a 
factor of Our data implies that this prediction is not 
valid at least for particles with St ^ 0.1. It remains to 
be checked whether the factor of \/2 difference between 
the radial and tangential rms relative speeds would be 
recovered at smaller Stokes numbers. At St *J 1, (w 2 ) 1 / 2 
and (w 2 ) 1 / 2 become exactly equal. 

The solid lines in Fig. [Jj correspond to the prediction 
of the PP10 model using eq. (|2"T1) for the angular aver- 
age of the trajectory structure tensor, Sri; . The equa- 
tion assumes that the direction of the particle separation 
R at any time is completely random, and predicts that 
(w 2 ) = (wp = i(w 2 ) for all particles (§3.2.2). This pre- 
diction is in good agreement with our simulation data. 
The equality of the radial and tangential rms speeds 
for St ^ 1 is expected, because the separation R of 
these particles at a friction timescale ago is significantly 
larger than the initial distance, r, and the direction of 
R is expected to be random with respect to r. On the 
other hand, the near equality of (w 2 ) 1 / 2 and {w 2 ) 1 / 2 at 
St ~ 0.1 — 0.2 is somewhat surprising. For r = I77, 
the backward separation of these particles does not con- 
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Fig. 7. — The rms relative speeds in the radial ((to 2 ) 1 / 2 ; cir- 
cles) and tangential ({ui 2 } 1 / 2 ; diamonds) directions. From top to 
bottom, data points show simulation results for particle pairs at 
r = 1, 0.5, and 0.25?;, respectively. Lin es a re predictions of the 
PP10 model. The solid lines adopt eq. f° r th e angular av- 

erage of Srij, which predicts that (m> 2 ) = (w 2 ) = ^(f 2 )- The 
particle separation behavior assumed here is exactly the same as 
in the solid lines in Fig. [6] for the 3D rms speed. The Richardson 
constant, g, is set to 1.6, 1.3 and 1.0 for r = 1, 0.5, and 0.25??, 
respec tivel y. The dashed and dotted lines for r = In are solutions 
of eq. (13 1 P for the radial and tangential speeds, respectively. The 
two lines reproduce the Saffman- Turner prediction for (i?? 2 ) 1 / 2 and 
{w^) 1 / 2 at small St. 

tribute to make the 3D amplitude, (w 2 ) 1 / 2 , of the relative 
velocity significantly larger than the S-T prediction. This 
suggests that the near equality of {w 2 ) 1 / 2 and (w 2 ) 1 ^ 2 is 
due to a conversion of the relative velocity from the tan- 
gential to the radial direction. The conversion is likely 
caused by the deviation of the particle trajectories from 
the fluid elements. Even though the deviation does not 
give rise to a considerable change for the 3D rms (w 2 ) 1 / 2 
at r = I77, it could efficiently alter the direction of w with 
respect to r. The trajectory deviation is stochastic, and 
thus tends to randomize the direction of w and equalize 
its radial and tangential components. This reduces the 
tangcntial-to-radial ratio below the S-T prediction. The 
randomization effect is expected to be more efficient in 
the slings events, where inertial particles are shot out of 
the streamlines of the flow, and encounter the trajecto- 
ries of other particles (Falkovich and Pumir 2007). At 
smaller r, the contribution from the backward separa- 
tion to the 3D relative speed is larger, and the random 
direction of the particle separation in the past also tends 
to equalize the radial and tangential relative speeds. 

When computing the solid lines, we used a bi- 
exponential form for $2, and the separation behavior 
adopted here is exactly the same as for the solid lines in 
Fig. [6] for the 3D rms speed. This means that the solid 
lines shown here correspond to those in Fig. [6] divided by 
\/3. The Richardson constant is set to 1.6, 1.3 and 1.0 
in the three lines for r = 1, 0.5, and 0.25?7, respectively. 

The dashed and dotted lines for r = lrj are the solu- 
tions of eq. (j3Tj) for the radial and tangential rms speeds, 
respectively. Eq. (|3"Tj) was derived from eq. (|26)) for 
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(5*Ty )ang, which assumes that the direction of the separa- 
tion change AR (rather than R itself) is random. When 
solving eq. (131[) . we used the same two-phase separation 
behavior (with g = 1.6) as in the corresponding solid 
line. At small St, the dashed and solid lines reproduce 
the S-T prediction that (w 2 ) — 2(w 2 ). The discrepancy 
between the simulation data and the S-T formula implies 
that the direction of R is more random than assumed in 
cq. flU for St in the range 0.1 £ St & 1. 

The dependence of the radial and tangential rms rela- 
tive speeds on r is similar to that of the 3D rms (see Fig. 
[6]). In an attempt to roughly fit them as power-law func- 
tions of r, we find that (w^) 1 / 2 and {w 2 ) 1 / 2 scale with r 
as oc r - 78 at St = 0.1 - 0.2. Similar to the 3D rms, the 
slight dips at St ~ 0.2 for both the radial and tangential 
relative speeds may be due to a numerical artifact. 

The simulations of Wang et al. (2000) found that 
the tangential-to-radial variance ratio, (w 2 )/(w 2 ), is ~ 
1.5 — 1.6 at St ~ 0.1 — 0.2. This is larger than the corre- 
sponding value (1.2 — 1.3) in our simulation, and closer 
to the S-T prediction. One possible reason for this dif- 
ference is that our flow has a much larger Reynolds num- 
ber. Although Wang et al. (2000) claimed that the ratio 
is independent of Re\ based on several simulations with 
Re,\ < 75, it is not clear if this is also true at Re\ 75. 
As speculated above, it is the deviation of the particle 
trajectories from the flow elements that tends to equalize 
the radial and tangential relative speeds of small parti- 
cles. Clearly, the trajectory deviation would be larger 
in flow regions with larger velocity gradients, where the 
flow experiences a faster velocity change. The proba- 
bility of finding large flow velocity gradients and hence 
large trajectory deviations increases with Re\. There- 
fore, the tangential-to-radial ratio is likely smaller at 
higher Reynolds numbers. Similarly, the frequency of 
sling events, which occur in regions with extreme flow ve- 
locity gradients, would increase with Re\ (Falkovich and 
Pumir 2007), and this also tends to reduce the tangential- 
to-radial ratio. However, we cannot rule out that the 
time resolution for the trajectory integration of the small- 
est particles in our simulation is not sufficient to allow 
an accurate measurement of (w 2 )/(w 2 ) at small St. 

6.1.3. Approaching and Separating Particle Pairs 

So far, our analysis for the particle relative velocity in- 
cluded all particle pairs at given distances. However, not 
all pairs at a small distance lead to collisions. Particles 
with a negative radial relative velocity, w v < 0, approach 
each other and may collide, while particles in pairs with 
w r > move away from each other. Since the goal of 
our study is to examine the collision of dust particles, it 
is appropriate to split the particle pairs at a given dis- 
tance into two groups: one with w v < and the other 
with w T > 0. We will refer to them as the minus and 
plus groups, respectively. Although only the first group 
is relevant for particle collisions, it is theoretically inter- 
esting to compare the relative velocity of particle pairs 
in the two groups. 

For the radial component, w T , of the relative velocity, 
we denote the variances in the minus and plus groups 
as and (w 2 ) + , respectively. In terms of the PDF, 

P(w T , St), of w T at a given St, the variances are written 

as (w r 2 )_ = J_ oo w 2 P(w r , St)dw r / J_ oQ P(w r ,St)dw r and 



( w 2) + = J°° w 2 P(w r , St)dw r / r o °° P(w z , St)dw r . We de- 
note the PDFs of the tangential component in the minus 
and plus groups as conditional PDFs, P(w t \w r < 0,St) 
and P(wt\w r > 0,St). The tangential variances in the 
two groups are then given by (w 2 ) T = J_ M w 2 P(wt\w r ^ 
0,St)dw t . Similarly, for the 3D amplitude, \w\, the 
minus and plus variances are expressed as (w 2 ) T = 
f™ oo w 2 P(\w\\w j: ^ 0,St)d\w\, where P(\w\\ Wj: ^ 0, St) 

are the PDFs of \w\ for approaching and separating pairs. 
In this subsection, we discuss the rms (or variances) of 
the particle relative velocity in the two groups, and the 
PDFs, P(w r ,St), P(w t \w T § 0,5*), and P(|u>|K § 
0, St), will be studied in §6.3. 

The data points in the left panel of Fig. ([8]) show the 
radial rms relative speeds of particle pairs at r = lr/ 
(squares) and 0.25?7 (circles). The right panel plots 
the tangential (squares) and 3D (circles) rms speeds at 
r = lr/. In both panels, the filled and open symbols cor- 
respond to particle pairs in the minus and plus groups, 
respectively, and the lines plot the overall rms relative 
velocities counting all particle pairs at given distances. 
If the velocity of St <C 1 particles closely follow the flow 
velocity, we expect that (w 2 )^ are determined by the 
variances of the longitudinal flow velocity increments, 
(Au 2 ) T , for the negative and positive Au T , respectively. 
The definition of (Au 2 )^ is given in Appendix B, and 
they correspond to the fluctuation amplitudes in the left 
and right wings of the PDF of Au r . In Appendix B, we 
find the ratio (Au 2 ) _ / (Au 2 ) + is 1.47 at the size of the 
computation cell, suggesting that, in the St — > limit, 

(w 2 )]/ 2 would be larger than (w 2 ) 1 ^ 2 by ~ 20%. The 
simulation result is consistent with this expectation. At 

St ~ 0.1, the ratio of (w 2 ) 1 / 2 to (w 2 ) 1 / 2 is ~ 25% for 
both r = 1 and 0.257? (see the left panel of Fig. |8]). Us- 
ing a similar analysis to the tangential component gives 

(w 2 )^ 2 ~ (Au 2 )^ 2 for St <C 1, where (Au 2 ) T are vari- 
ances of the transverse flow velocity increment, Au t , con- 
ditioned on the sign of longitudinal increment Au r (see 

Appendix B). The ratio (w 2 )_/(w 2 )V 2 is found to be 
1.16 at St = 0.1—0.2, consistent with the ratio of (Au 2 )^ 
to (Au 2 )+ at the cell size (see Appendix B). For the the 

3D amplitude ((w 2 ) 1 ^ 2 ), the rms ratio between the minus 
and plus groups is 1.18 at St <C 1. 

As St increases, the relative speed for the plus group 
first decreases slightly and reaches a minimum at St ~ 
0.4 in all cases with r = \r\. This can be explained by 
considering the effects of the particle memory and the 
particle separation backward in time. Particle pairs in 
the plus group with w r > are coming from smaller dis- 
tances, meaning that the separation of the particles was 
smaller in the near past. As St increases from 0.1 to 0.4, 
the contribution from the particle memory of the flow ve- 
locity difference becomes relatively more important, and 
this contribution tends to reduce the relative speed since 
the particle distance was smaller in the immediate past. 
However, if we look back further into the past (i.e., at 
larger \t\), the two particles may pass each other, and 
their distance would make a transition from decreasing 
to increasing as |r| keeps increasing. This explains the 
increase of (w 2 ) + , (w 2 } + , and (w 2 ) + at St <; 0.4. The 
minimum of (w 2 }+ for r = 0.2577 appears at St = 0.2 
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Fig. 8. — The rms relative speeds for particle pairs approaching (w r < 0; filled symbols) or separating (ui r > 0; open symbols) from 
each other. Lines correspond to the overall rms relative speeds counting all particle pairs. Left panel: The radial rms speed (w^)]^ 2 ) in 
the minus and plus groups with r = lr] (squares) and 0.25r; (circles). Right panels: the tangential (squares) and 3D (circles) rms relative 
speeds for particle pairs at r = lr). 



instead of St — 0.4. This is because, for smaller r, it 
takes a shorter time for the particle distance in the past 
to change from decreasing to increasing. 

For approaching particles in the minus group, the par- 
ticle distance would increase monotonically toward the 
past. Therefore, the relative speed for the minus group is 
expected to monotonically increase with St for small par- 
ticles with r p <, Tl. This is confirmed by the filled data 
points in Fig- [HI except the slight dips at St ~ 0.2. These 
slight dips are not expected, and again may be caused by 
insufficient numerical accuracy in the trajectory integra- 
tion for the smallest particles in our simulation (§6.1.1). 
Fig. [5] shows that approaching particles tend to have a 
larger relative speed than the separating pairs. The dif- 
ference between the two groups first increases with St 
and then starts to decrease at St > 1. At St ^ 6.2, the 
rms relative speeds in the two groups are close and equal 
to the overall rms. The reason is that, for these larger 
particles, the particle separation at a friction timescale 
ago becomes independent of the "initial" condition at 
t = 0. 

The rms relative velocity in the minus group is found to 
be larger than the overall rms for particles St ^ 6, and 
an immediate implication is that using the overall rms 
to calculate the collision energy may lead to an under- 
estimate for these particles. This has not been pointed 
out by previous studies, which considered the overall rms 

only. The difference between the 3D rms, (w 2 ) 1 / 2 , in 
the minus group and the overall rms {w 2 ) 1 / 2 peaks at 
St = 1, where (w 2 )]/ 2 is larger by 25%. Therefore, the 
collision energy may be underestimated by up to 60% 
if one uses the overall rms relative speed. The predic- 
tion of the PP10 model was for the overall rms, and it 
could be modified to predict the relative velocity specif- 
ically for approaching particle pairs. For that purpose, 
the backward separation behavior of approaching parti- 
cle pairs must be specified and taken into account. We 

also find that (w 2 ) 1 / 2 and (w 2 ) 1 / 2 almost coincide for 
all St, meaning each component of the relative velocity 



provides equal amount of collision energy. On average, 
the radial component provides 1 /3 collision energy, while 
the remaining 2/3 is contributed by the two tangential 
components. 

6.2. The Collision Kernel 

The prediction of the collision rate is one of the main 
goals of our study of the particle relative velocity. If the 
mean number density of inertial particles of a given size 
is rip, then the collision rate per unit volume between 
these identical particles is estimated by \n^T : where T 
is the collisional kernel and the factor 1/2 is used to 
avoid counting the same pair twice. The collision ker- 
nel depends on the particle cross section and the particle 
relative speed. The collision rate is also affected by the 
spatial clustering of inertial particles in turbulent flows. 
Due to their finite inertia, inertial particles do not ex- 
actly follow the flow velocity and have been found to ex- 
hibit inhomogeneous distribution even in incompressible 
turbulence. Turbulent clustering of inertial particles has 
been extensively investigated in the literature (see e.g., 
Maxey 1987, Sundaram and Collins 1997, Wang et al. 
2000, Cuzzi et al. 2001, Hogan & Cuzzi 2001, Balkovsky 
et al. 2001, Zaichik et al. 2003, Falkovich and Pumir 2004, 
Cuzzi et al. 2008, Pan et al. 2011). The general physical 
interpretation for turbulent clustering is that the vorti- 
cal structures in turbulent flows tend to expel inertial 
particles. The particles are pushed out of high-vorticity 
regions by the centrifugal force, leading to the formation 
of clusters in strain-dominated regions. The degree of 
turbulent clustering is usually quantified by the so-called 
radial distribution function (RDF), g(St,r). The RDF 
is defined such that the number of particles in an in- 
finitesimal volume, dV, located at a distance r from a 
given particle, is n p g(St, r)dV. The RDF is essentially 
the enhancement factor for the probability of finding a 
neighboring particle. The contribution of turbulent clus- 
tering to the collision kernel (or the collision rate) for 
identical particles is a factor of g(St, d p ) at a distance of 
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Fig. 9. — The radial distribution function, g(St, r), as a function 
of the Stokes number at r = 1 (solid), 0.5 (dashed) and 0.25r? 
(dotted). 



the particle diameter, d p . 

Fig. [9] plots the RDF as a function of St at three dis- 
tances, r = 1, 0.5 and 0.25??. Consistent with previous 
studies, the RDF is largest for St ~ 1 particles, whose 
friction timescale couples with the smallest scale of the 
turbulent flow. For particle distance, r, in the dissipation 
rage, the RDF increases toward smaller r as a power-law, 
i.e., g(St,r) oc . The scaling exponent /j peaks at 
St 1, and approaches zero in the limits St « 1 and 
St 3> 1. We measured /x using the values of r shown in 
Fig. ij and found that fx = 0.73 for St = 0.78, which 
is consistent with the result of Pan et al. (2011). The 
interested reader is referred to Pan et al. (2011) for the 
scaling exponent, /i, as a function of St. Existing coag- 
ulation models for the growth of dust particles in pro- 
toplanetary disks usually ignore the effect of turbulent 
clustering. This would underestimate the collision rate, 
especially for dust particles with St ~ 1. At r = 0.25??, 
the RDF is already ~ 20 for St ~ 1. The size, a p , of dust 
particles is typically much smaller than the Kolmorgorov 
length scale, r\ (~ 1 km), of the protoplanetary turbu- 
lence. Due to the large separation between r\ and o p and 
the power-law increase of the RDF toward small r, the 
effect of turbulent clustering would provide a significant 
contribution to the collision rate. Before discussing the 
clustering effect in more details, we will first consider the 
formulation of the collision kernel and the contribution 
from turbulence-induced collision speed. 

Saffman and Turner (1956) presented two formulations 
for the collision kernel. The formulations were based 
on spherical and cylindrical geometries, respectively, and 
were thus named as the spherical and cylindrical formula- 
tions by Wang et al. (1998). In the spherical formulation, 
the collision kernel is written as T sph — Anrd^g{St, d p )F~ . 
Here F~ represents the radial inward flux from a dis- 
tance of dp to a given particle. In terms of the PDF, 
P(w T ,St), of the radial relative speed, the inward flux 

is given by F~ = — f_ w x P(w T , St)dw T , In statisti- 



cally isotropic turbulence, as is the case for our simu- 
lated flow, the particle statistics is expected to become 
isotropic when the dynamics is fully relaxed. In that 
case, the inward flux, F~ , is equal to the outward flux 
F + (= w T P(w r , St)dw T ), which is confirmed by our 
simulation data. This is because the average radial ve- 
locity (w T ) = from isotropy and (u> r ) = F+ — F~ by 
definition. We thus have F r + = F r ~ = ^{\w T \) where 
(|u>r|) is the ensemble average of the absolute value of 
the radial relative velocity. The collision kernel can then 
written as T sph = 2nd%g{\w I \). 

In the left panel of Fig. (fTD|) . we plot the simulation 
result for (\w T \) at three particle distances r = 1,0.5, 
and 0.25?7. For comparison, we also show the data (cir- 
cles) for the rms of the radial relative speed at r = lrj. 
Qualitatively, (|u> r |) as a function of St and r is similar 
to the rms. It is smaller than the rms, as it corresponds 
to the lst-order moment of the PDF, P(w Y ,St), of w r . 
Most theoretical models, including our own, for the par- 
ticle relative velocity are based on the computation of 
the variance (e.g., (w 2 )) of the relative speed, and can- 
not be directly applied to predict (|iw r |). The conversion 
between (w 2 ) and (|u> r |) relies on the PDF of w T , which is 
difficult to predict. We thus did not attempt to fit (\w T \) 
with a model prediction. The simulation result for the 
PDF, P(w Y ,St), as a function of St will be discussed in 
§6.3. 

Similar to the S-T formula (eq. ITT]) for the variances 
of the relative velocity, we would predict that (\w T \) — 
(|Au r |) in the limit St -C 1, where (|Ait r |) is the abso- 
lute average of the longitudinal flow velocity increment. 
At i ^ 77, (|Au r |) is expected to scale linearly with I. 
We find that, for small particles with St ^ 0.4, (|tu r |) 
at r = I77 is consistent with (|Au r |) at £ = lrj, which is 
~ 0.19-1/,,) in the simulated flow. At smaller r, { | w r | ) is 
found to scale with r as oc r 9 , which is slightly shal- 
lower than the linear scaling. This is likely caused by 
the contribution from the effect of slings or caustics at 
small r. The scaling is steeper than r 78 for the radial 
rms velocity (see §6.1.2), suggesting that (|u> r |) reflects 
the flow velocity scaling better. The behavior of (|Au r |) 
at small St also appears to be more regular than that 
of {w^) 1 / 2 . For all three values of r, the (|w r |) curves 
become flat at St < 0.4, as expected from the prediction 
(|u> r |) = (|Au r |). The likely reason is that (|ie r |) rep- 
resents statistics at a lower order than the rms, and is 
thus less affected by the rare and extreme sling events, 
or by the numerical uncertainty in the particle trajectory 
computation. 

The inset of the left panel shows the ratio of the abso- 
lute average to the rms. The ratio depends on the shape 
of the PDF of u> r , and particularly on the central part 
of the PDF because both (|u> r |) and the rms are lower- 
order moments. As a reference, if the PDF P(w r , St) is 
Gaussian, we have (K|)/(u£) 1/2 = (2/tt) 1 / 2 = 0.80, and 
for an exponential PDF it is equal to l/y/2. The central 
part of an exponential PDF is sharper than the Gaussian 
case. Generally, the ratio decreases as the PDF becomes 
fatter. As seen in the inset, the ratio changes with both 
St and r. For r = lrj, the ratio is ~ 0.68 at St ~ 0.1-0.2. 
With increasing St, it decreases and reaches a minimum 
of 0.45 at St ~ 1, indicating a highly non-Gaussian PDF. 



22 



3 




10- 1 10° 10 1 10 2 
St 



10- 1 



3D rms 



10" 



10 u 



10 1 
St 



10 2 



10 u 



G 


o o 










1.0 








: 


0.8 






0.6 




r - 


0.4 






0.2 
10 






- 1 10° 


10 1 10 2 
St 


10 1 

St 




10 2 



Fig. 10. — Left panel: The average of the absolute value of the radial relative speed, (|«i r |), as a function of the Stokes number. Solid, 
dashed and dotted lines correspond to the particle distance r = 1,0.5 and 0.25ri, respectively. For comparison, the circles show the rms, 
(id 2 ) 1 / 2 , of the radial relative speed at r = lr). The inset plots the ratio, (\w T \) / (w^) 1 ^ 2 , at r = 1 (solid), 0.5 (dashed) and 0.25?; (dotted). 
Right panel: same as the left panel, but for the 3D amplitude, \w\, and the 3D rms, (u; 2 ) 1 ' 2 , of the relative velocity. 



As will be shown in §6.3, as St increases to ~ 1, the fat- 
ness of the PDF P(w T , St) increases: its central part be- 
comes sharper while the tails get broader. As St further 
increases above 1, the trend for the fatness of the PDF 
shape is reversed, and the ratio (\w?\) / (w 2 ) 1 / 2 increases. 
At St ~ 800, it reaches 0.78, close to the expected value 
for a Gaussian PDF. The ratio also decreases with de- 
creasing r for St < 6. This is due to the sharpening 
of the central part of the PDF with decreasing r (see 
§6.3.3). At larger St, the PDF shape and hence the ratio 
{\w T \) / {w 2 ) 1 ^ 2 are independent of r (§6.3.3). 

We next discuss the cylindrical formulation, in which 
the collision kernel is written as r cyl = ird 2 g(\w\), with 
(\w\) being the average of the 3D amplitude of the rel- 
ative velocity. This formula assumes that all particles 
inside a cylinder of length (\w\)dt located at a distance 
dp from a given particle will collide with the particle in a 
time interval dt. Our simulation result for (\w\) is shown 
in the right panel of Fig. (fTU|). which is very similar to 
the left panel for { | w r | ) . At small St <, 0.4, we find (\w\) 
also scales as r 09 with r. The ratio of (\w\) to the rms, 
(iu 2 ) 1 / 2 , also shows a dip at St ~ 1 and approaches 0.9 
at the largest St. Similarly, this can be understood from 
the PDF of | «> | as a function of St (see §6.3). The value 
of 0.9 is expected from a 3D Gaussian distribution. 

In Fig. [11] we plot the products of the RDF and the ab- 
solute average of the relative velocity, 2g(St,r)(\w I \) /u n 
(solid lines) and g(St, r)(\w\)/u v (dashed lines), for the 
spherical and cylindrical formulations of the collision ker- 
nel, respectively. The two products correspond to r sph 
and r cy measured at a distance of r and normalized to 
irdpU^. We refer to the products as the normalized kernel 
or the collision kernel per unit cross section. The lines 
from bottom to top are measured at r = 0.25, 0.5 and lr], 
respectively. At each r, the solid and dashed lines almost 
coincide, meaning that r sph and T cyl are nearly equal at 
all St and r. This suggests that (\w T \) ~ 0.5(|io|) since 
pprypcyl = 2{\w T \)/(\w\). The two collision kernels have 
a noticeable difference only at St = 0.1 — 0.2, where r sph 



is smaller than r cyl by <^ 5%. This is consistent with 
the result of Wang et al. (2000), who found that the pre- 
dictions of the two formulations are equal except a differ- 
ence of £ 25% in the limit St -> 0. Wang et al. (2000) 
also showed that the spherical formulation provides an 
almost exact description for the particle collision rate. 
The near equality of r sph and r cyl at all St and r sug- 
gests that one can apply either formulation to evaluate 
the collision rate. 

Fig.rn]shows that 2g(St,r){\w T \) and g(St,r)(\w\) are 
independent of r for St > 1. From Figs. [5] and [T0l 
we see that, with decreasing r, the RDF, g(St,r), in- 
creases, while the absolute average of the relative ve- 
locity decreases. Interestingly, for St > 1, the increase 
of g(St, r) almost exactly compensates the decrease of 
(\w T \) (or (|it>|)). For example, at St = 1.55, we have 
g{St,r) oc r ~ - 55 , while both (\w T \) and (\w\) scale with 
r asoc r 57 . As a result, their products are independent 
of r. According to the model of Wilkinson et al. (2006) 
and Gustavsson & Mehlig (2011), the effect of caustics 
provides a scale-independent contribution to the normal- 
ized collision kernel. At a given St, the caustic contribu- 
tion would become dominant at sufficiently small r. The 
r-independence of the normalized kernel for St > 1 parti- 
cles in our simulation suggests that the effect of caustics 
for these particles is already dominant at r <J 77. This in- 
dependence also implies that, for St > 1, one can safely 
apply our result for 2g(St,r)(\w T \) (or g(St, r)(\w\)) at 
r ~ X] to estimate the collision kernel at the particle di- 
ameter, dp, even though, in the case of dust grains in 
protoplanetary disks, d p is smaller than 77 by orders of 
magnitude. 

We find that the collision kernel per unit cross section 
is almost constant for St between 1 and the peak Stokes 
number, St m ~ 30. This is because of the opposite trends 
of the RDF and the relative velocity with St in this range. 
The products 2g(St,r)(\w T \) and g(St, r)(\w\) increase 
only by 50% as St goes from 1 to 30. The kernel starts to 
decrease at St ~ 30, and finally scales with St as St~ x / 2 
in the large particle limit with r p ^> Tl. Interestingly, 
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Fig. 11. — The collision kernels normalized to 7rdpU, ; . The solid 
lines plot 2g(St, r) (\w T \) / for the spherical formulation, and the 
dashed lines show g(St, r){|io|)/utj for the cylindrical formulation. 
The lines from bottom to top correspond to normalized kernels 
measured at r = 0.25, 0.5 and I77, respectively. 

in the small particle limt the normalized kernels experi- 
ence a shape rise as St increase from 0.1 to 1 (see, e.g., 
Svmdaram & Collins 1997), suggesting that the collision 
frequency greatly accelerates once the particle size grows 
past St ~ 0.2. For r J> 0.25?y, the rise is mainly con- 
tributed by the fast increase of the RDF toward St ~ 1. 
Based on the model of Wilkinson et al. (2006), the abrupt 
rise of the normalized kernel corresponds to the rapid 
formation of caustics, and it was claimed that the effect 
of caustic formation may be modeled as an activation 
process. If we fit the normalized kernel by exp(— A/ St) 
(Wilkinson et al. 2006, Falkovich & Pumir 2007), the 
activation value, A, is found to be around 0.8. 

We give a simple physical explanation for why the sling 
effect or the caustic formation causes the rapid rise of 
the normalized kernel at St ~ 1. The sling events occur 
at places where the flow velocity gradient is larger than 
r" 1 (Falkovich et al. 2002). At these locations, the flow 
velocity changes at a timescale faster than the response 
of the particle. For small particles with r p <C r ?) , the 
probability of the sling events corresponds to the high 
tails of the probability distribution of the flow velocity 
gradient. With increasing r p , the probability becomes 
larger as it samples toward the inner parts of the flow 
gradient PDF. Because the rms flow gradient is on the 
order of t" 1 , the probability experiences a fast increases 
as Tp" 1 approaches r" 1 , leading to a sharp rise in the 
sling frequency and hence in the collision rate at St ~ 1. 

At St < 0.8, the increase of the RDF toward small r 
is slower than the decrease of the relative speed. Conse- 
quently, the normalized kernel decreases with decreasing 
r, as seen in Fig. 1111 This r-dependence needs to be con- 
sidered when evaluating the collision kernel at d p <C 
for St < 0.8 particles. A simple approach is to mea- 
sure the scaling exponent of the normalized kernel with 
r and then extrapolate the kernel to d p . For example, 
the normalized kernel scales with r roughly as r 68 at 



St = 0.1. However, the validity of this method would fi- 
nally break down as r keeps decreasing. At smaller r, the 
contribution of caustics becomes more important. The 
theoretical result of Gustavsson & Mehlig et al. (2011) 
indicates that, at a given St, there exists a critical parti- 
cle distance, r c , below which the effect of caustics dom- 
inates. Since the caustic contribution is r-independent, 
the normalized kernel would become constant at r < r c . 
The critical distance r c is generally a function of St. Our 
simulation result for the normalized kernel suggests that 
r c 1^7 for St > 0.8 particles (see above), and r c <, 0.25r? 
for St < 0.8. If the dependence of r c on St and the 
scaling behavior of the normalized kernel above r c were 
determined, one could estimate the collision kernel of 
St < 0.8 particles for any value of d p . However, a numer- 
ical study to investigate r c and the caustic contribution 
as a function of St, for St < 0.8 particles, is computation- 
ally challenging, as it requires a huge number of particles 
to measure statistics of the RDF and the relative speed 
at r -C rj. A possible compromise would be to simulate 
the turbulent flow at lower resolution in order to allow 
the computation of a much larger number of particles. 
We defer such a study to a future work. 

In existing coagulation models of dust particle growth 
in protoplanetary disks, the collision kernel is typically 
set to r coa = Trd^w 2 ) 1 / 2 (e.g., Dullemond and Dominik 
2005), which is of the cylindrical type, but different from 
F cyl defined above. It ignores the effect of turbulent clus- 
tering and uses the 3D rms instead of the average 3D 
amplitude, (|u>|), of the relative velocity. Clearly, r cyl is 
larger than r coa by a factor of g(St, r p )(\w\) /(w 2 ) 1/2 . As 
seen in Figs .191 and ITUl g(St,r p ) is larger than 1, while the 
second factor (w 2 ) 1 ^ 2 / (\w\) is smaller than 1. From our 
simulation data, we find that g(St,r p )(\w\) /{w 2 ) 1 ^ 2 is 
close to unity for St <; 6 or St ^ 0.2, and the difference 
between r coa and r cyl in this St range is within 50%. 
The agreement of r coa and T cyl at St -C 1 is simply a 
coincidence, where g(St,r p ) happens to compensate the 
ratio {\w\) / {w 2 ) 1 ' 2 . On the other hand, for St between 
0.2 and 6, r cyl is significantly larger than r coa due to 
the large values of g(St,r p ). The difference is largest 
at St ~ 1. For example, at St = 0.78, r coa is smaller 
than r cyl by a factor of 3.8 and 6.9 for r = lrj and 0.25??, 
respectively. The difference would be even larger as r fur- 
ther decreases. We thus conclude that the collision kernel 
commonly used in dust coagulation models applies only 
for St <; 6 or in the limit St — > 0, and it significantly 
underestimates the collision rate at St ~ 1 . 

Besides the imprecise formulation of the collision ker- 
nel, coagulation models for dust particles in protoplane- 
tary disks also adopt estimates of the rms collision veloc- 
ity induced by turbulence that are inaccurate, especially 
for particles with St <, 6. The range St <, 6 may include 
all particle sizes up to chondrulcs for typical disk param- 
eters at a radius of 1 AU. The commonly-adopted rms 
collision velocities are based on the model of Volk et al. 
(1980) (or its later developments, e.g., Ormel and Cuzzi 
2007). There are two main problems in using these mod- 
els to estimate the collision rate. First, as mentioned ear- 
lier, the Volk et al. model overestimates the peak relative 
speed at St m (§6.1), which would result in overestimating 
the collision rate by a factor of two or so for inertial-range 
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particles and/or large particles with St > St m . Second, 
the model of Volk et al. does not reflect the dependence 
of the relative speed on the particle distance, r, because 
it does not keep track of the particle separation (PP10). 
Therefore, the model does not account for the S-T regime 
and may be viewed as corresponding to r — > 0. The pre- 
dicted relative speed by Volk et al.'s model sharply drops 
to zero at St = 1, suggesting the collision kernel is ex- 
actly zero at St < 1. This is clearly not the case from 
the simulation result. The relative velocity does have a 
dependence on r for St <> 6, which, together with the 
r-dependence of the RDF, is crucial for determining the 
collisional kernel for St <, 6 particles. The cancellation of 
the two dependences gives a finite collision kernel around 
St ~ 1. These features cannot be captured by the Volk 
et al. model, which is thus likely inapplicable to the es- 
timate of the collision rate for St <, 6 particles. 

6.3. The PDF of the Particle Relative Velocity 

In this section, we explore the probability distribution 
of the relative velocity of nearby inertial particles. An 
accurate estimate of the PDF of the particle relative ve- 
locity is important for modeling the growth and evolution 
of dust particles in protoplanetary disks. As mentioned 
in the Introduction, the outcome of particle collisions de- 
pends on the collision velocity, and due to the random na- 
ture of the turbulent-induced relative velocity, collisions 
of particles with exactly the same properties may have 
different outcomes, and thus using a single value, e.g., the 
rms, for the collision speed of particles of a given size is 
insufficient. The probability distribution of the collision 
velocity is needed to calculate the fractions of collisions 
resulting in sticking, bouncing or fragmentation. 

The physical picture of PP10 shows that the relative 
velocity of inertial particles depends on the flow veloc- 
ity difference the particle pair saw within a memory 
timescale or so. This suggests that the statistics of the 
velocity difference in the carrier flow is crucial for the 
understanding of the relative velocity PDF of inertial 
particles. Therefore, we analyzed the PDFs, P u (Au T ,£) 
and P u (Au t , I), of the longitudinal and transverse veloc- 
ity increments, Au r and Au t , as functions of the length 
scale, £, in our simulated flow. The results are discussed 
in details in Appendix B. The flow velocity PDFs are 
used in the physical explanation for the PDF of the par- 
ticle relative velocity as a function of particle inertia in 
§6.3.2. 

6.3.1. The PDFs of the radial and tangential components 

In Fig. [HI we show the PDFs, P(w T ,St), of the ra- 
dial component of the relative velocity as a function of 
the Stokes number. All the PDFs are measured from 
particle pairs at a distance of lr/. The relative speed is 
normalized to the Kolmorgorov velocity, u n , and each 
PDF is normalized to its value at the central peak. In 
the left panel, the thin dashed line corresponds to tracer 
particles (St — 0). The shape of this line is found to 
be identical to the PDF, P u (Au T ,£), of the longitudi- 
nal flow velocity increment, Au r , at the computational 
cell size (£ — 1.77/; see the top line in the left panel of 
Fig. [TH]in Appendix F3). This is expected as the tracer 
particles exactly follow the flow velocity, and P u (Au T ,£) 
is independent of I in the dissipation range. The solid 



color lines of increasing width show the PDFs of particles 
with increasing St. This corresponds to the increase of 
the rms relative speed with St for St Ss 24.9 (see Figs. 
[7] and©. For St < 1.55, the PDF of w T has a negative 
skewness, which is inherited from the flow velocity PDF 
P u (Au r ,£). The PDF becomes symmetric at St 3.21. 
It is interesting to note that, as St increases, the tails 
of the PDFs become broader, while the innermost part 
remains unaffected and equal to the PDF of the tracer 
particles. As to be shown in §6.3.2, the amplification of 
the PDF tails for small particles correspond to the ef- 
fect of slings. Due to the tail amplification, the overall 
PDF shape becomes fatter as St increases to ~ 1.55. E3 
With increasing St, the amplification effect proceeds to- 
wards the inner parts of the PDF, leading to a sharp 
cusp-like shape at the center, especially for St > 3.11. 
For particles with St J> 3.11, we see that the slope of the 
outer parts of the PDF tends to steepen when extending 
to higher tails, i.e., the PDF shape becomes thinner at 
larger values of w v . This thinning trend toward the high 
tails causes a decrease in the overall fatness of the PDF 
for St above 3.11. 

In the right panel, the PDF becomes narrower as the 
Stokes number increases above 49.7, corresponding to the 
decrease of the rms velocity with St in the large particle 
limit (Figs. [7] and [5]). The friction timescale of these par- 
ticles is larger than the correlation timescale (XL) of the 
flow velocity at the largest scales, meaning that the mem- 
ory time of the flow velocity is shorter than the memory 
of the particles. This induces a factor of Xl/t p in the 
variance of the relative velocity (§3.2.4), causing a de- 
crease of the PDF width at large St. The dotted line 
in the right panel is the Gaussian fit to the PDF of the 
largest particles (St = 795) in our simulation. The fric- 
tion timescale of these particles is 54 times larger than 
Xl, suggesting that the assumption of a Gaussian relative 
velocity PDF applies only in the extreme limit r p ^> Xl. 

Fig. EH shows the PDF, P(w t \w r < 0,t p ), of a tan- 
gential component of the relative velocity conditioned 
on w r < 0. The measurement of P(wt\w T < 0,t p ) only 
counts particle pairs approaching each other. The figure 
is plot in the same way as Fig. [T2 for the radial com- 
ponent. Again, the thin dashed line in the left panel 
is for tracer particles (St = 0). It corresponds to the 
PDF of the transverse difference, Aut, of the flow veloc- 
ity conditioned on Au r < 0. Our data confirms that the 
shape of P(wt\w r < 0,0) for tracer particles at r = lr) 
is close to the dashed line in the right panel of Fig. [IS] 
for P u (Au t \Au r < 0,£) at I = l.7r) (Appendix B). The 
qualitative behavior of P(wt\w I < 0,St) as a functions 
of St is similar to that of P(w T , St). 

For St = 0.19 particles, the tails of the radial PDF are 
significantly amplified with respect to tracers (Fig. [T2|) , 
while the conditional PDF of the tangential component 
almost coincides with the tracer PDF (Fig. fTS")) . From the 
physical picture for the PDF behavior given in §6.3.2, the 
effects of the particle memory and the backward separa- 
tion tend to amplify the PDF tails of St = 0.19 particles. 
However, for the tangential PDF, this amplification effect 

11 For definiteness, throughout the paper we use "fat" or "thin" 
to describe the shape of the PDF. The fatness can be quantified, 
e.g., by kurtosis. On the other hand, the extension or width of 
the PDF, corresponding to the rms, is described as "broad" or 
"narrow" . 
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Fig. 12. — The PDF of the radial component (w T ) of the relative velocity at r = Irj as a function of the Stokes number. The relative 
speed is normalized by the Kolmorgorov velocity, u v , and each PDF is normalized to its peak value at w T = 0. The left panel shows the 
PDFs for particles with St < 24.9, while the right panel shows results for larger particles with St > 49.7. The dotted line in the left panel 
is the PDF of the radial relative speed of tracer particles (i.e., St = 0). The dashed line in this panel is the stretched exponential function 
with a = 4/3, which provides a good fit for the PDF tails of St = 24.9 particles. In the right panel, the dashed line corresponds to the 
Gaussian fit to the largest particles in our simulation. 
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Fig. 13. — The PDFs of the tangential component, Wt, of the relative velocity conditioned on w r < 0. The PDFs are measured from 
particle pairs approaching each other. The normalization of the PDFs is the same as in Fig. 1121 The left and right panels show results 
for particles with St < 24.9 and St > 49.7, respectively. The dotted line in the left panel is the PDF of the tangential relative velocity of 
tracer particles (St = 0) conditioned on w T < 0. The dashed line in this panel is the stretched exponential function with a = 4/3. The 
dashed line in the right panel corresponds to the Gaussian fit for St = 795 particles. 



is counteracted by the conversion of the relative veloc- 
ity from the tangential direction to the radial direction, 
which reduces the PDF width of the tangential compo- 
nent. As discussed in §6.1.2, the conversion is caused by 
the deviation of inertial particle trajectories from flow 
elements and the particle memory of the flow velocity 
in the past, which randomize the direction of the rela- 
tive velocity relative to the particle separation, r. The 
conversion is expected to especially efficient at the PDF 
tails, corresponding to the sling events. It appears that 
the two opposite effects cancel out for the tangential PDF 
of St ~ 0.19 particles, as it almost coincides with the 
dotted line for tracers. On the other hand, both effects 



broaden the PDF of the radial component, leading to 
significantly amplified tails with respect to tracers. 

Unlike P(w T ,St), which has a negative skewness for 
St < 1.55 particles, P(wt\w T < 0,St) is symmetric at 
all St. As mentioned earlier, the symmetry of the PDF 
in a tangential direction is expected from the statistical 
isotropy. We find that the left wing of the radial PDF 
P(iVr, St) almost coincides with that of P(wt\w r < 0, St) 
at all St. This is apparently due to the randomization 
of the relative velocity direction discussed above. On the 
other hand, the right wing of P(w r , St) is narrower than 
that of P(w r \w r < 0,St) until it becomes symmetric at 
St > 3.11. For particle collisions, we are mainly inter- 
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ested in the PDFs for approaching pairs, i.e., the left 
wing of P(w T , St) and the entire PDF, P(w t \w T < 0, St). 

We give a more quantitative description for the shape 
of P(wt\w r < 0, St). The description also applies to the 
left wing of P(w T , St), as it coincides with the left wing 
of P(wt\w T < 0,St). We first quantify the fatness of 
the tangential PDF by computing the kurtosis, defined 
as (wf) - 1 (w 2 ) 2 _ , for P(w t \w r < 0,t p ). At r = lr), the 
kurtosis for St = 0.1 and 0.19 particles is ~ 11, which is 
already much larger than 3 for a Gaussian PDF. With 
increasing St, the kurtosis first increases due to the tail 
amplification. It reaches a maximum value of 36 at St = 
0.78, indicating extremely high non-Gaussianity. The 
kurtosis decreases slightly to 32 at St — 1.55, and then 
drops rapidly and approaches ~ 3 at the largest St(= 
795) in our simulation. This decrease corresponds to 
the thinning trend of the high PDF tails for these large 
particles. We also measured the kurtosis for the PDFs 
at smaller r, and found that it increases with decreasing 
r for St 6.2. The PDFs are fatter at smaller r because 
the effect of the tail amplification on the PDF shape is 
relatively stronger (see discussions in §6.3.3). Also for a 
smaller r, the maximum kurtosis occurs at smaller St. 

Following Sundaram & Coitions (1997) and Wang 
et al. (2000), we attempted to fit P(w t \w I < 0,r p ) 
with stretched exponential distributions. The generic 
stretched exponential function is defined as, 
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where T is the Gamma function. The variance of P sc is 
given by /3 2 T(3/a)/T(l/a). Thus, to fit a given PDF by 
eq. (f34|) with a chosen value of a, (3 can be fixed by the 
variance of the PDF. The index a controls the shape of 
the PDF, and smaller a corresponds to fatter tails. The 
PDFs for St = 0.1 and 0.19 at r = lrj almost have the 
same shape, and both can be fit by a stretched exponen- 
tial function with a = 0.67. This value of a is consistent 
with that (0.7) used to fit the PDF tails of the flow ve- 
locity difference at £ = 1.7 (see Appendix B and Fig.[T8|). 
At St = 0.39, 0.78 and 1.55, the best-fit values of a are 
0.52, 0.48, and 0.49, respectively. The decrease of a in 
the St range from 0.1 to 0.78 also indicates increasing 
fatness of the PDF. The PDF shape at St = 1.55 is very 
close to that for St = 0.78 particles. 

For particles with 3.11 < St < 49.7, the PDF cannot 
be well fit by a single stretched exponential function. The 
PDFs of these particles are more complicated, due to the 
existence of sharp cusps at the center and the steepening 
trend of the PDF slope toward to the far tails. These 
features cannot be captured simultaneously by a single 
stretched exponential. It is, however, possible to fit these 
PDFs with a combination of two different stretched func- 
tions for the cusp and the tails respectively. We postpone 
a detailed study of the fitting functions for these inter- 
mediate particles to a future paper with simulations at 
a higher resolution. To give a quantitative idea for the 
PDF shape of these particles, here we simply list the best- 
fit a for the far tails without accounting for the central 
cusp. At St = 3.11, 6.21, 12.4, 24.9 and 49.7, a for the 
PDF tails is found to be 1, 1.1, 1.3, 1.33, and 1.45, re- 
spectively. The stretched exponential fits for the PDF 
tails of St = 24.9 particles are shown as dotted lines in 



Fig. [TS] and [T3J where the index a is set to 4/3. Such 
a 4/3 stretched exponential PDF was predicted by Gus- 
tavsson et al. (2008) assuming an exact Gaussian flow 
velocity field with Kolmorgorov scaling and a rapid tem- 
poral decorrelation. An alternative derivation for the 4/3 
stretched exponential is given in §6.3.2 using the physi- 
cal picture of the PP10 model. Our derivation does not 
assume a short temporal correlation for the flow velocity, 
and is thus more general than that of Gustavsson et al. 
(2008). 

Starting from St — 99, the central cusp becomes suf- 
ficiently small, leading to simpler PDF shapes. This al- 
lows the entire PDF to be satisfactorily fit by a single 
stretched exponential again. The measured a values for 
St = 99, 199, 397 and 795 are 1.5, 1.65, 1.75 and 1.9, 
respectively. Note that the PDF for St — 795 is close 
to Gaussian, but the best-fit value for a is actually 1.9 
instead of 2. A similar trend of the PDF fatness and 
the best-fit a as a function of St was found in previ- 
ous studies with low-resolution simulations (Sundaram 
& Collions 1997; Wang et al. 2000) 

6.3.2. Physical picture for the PDF behavior 

We give a detailed explanation for the behavior of the 
relative velocity PDF based on the physical picture of the 
PP10 model. The picture was illustrated in Fig. [TJ We 
first consider particles with r p ^ TL. In §3.2.4, we showed 
that the 3D variance of the relative velocity of these 
particles may be roughly estimated by (w 2 ) ~ Su (Rp), 
where S!y is the flow structure tensor and R p is the pri- 
mary distance. For simplicity, we have neglected the ef- 
fect of the temporal correlation function, <I>2, which may 
provide a factor (T(R p ) /r p ) of order of unity for particles 
with t p <, Tl. R p was estimated by R 2 = r 2 + (w 2 )r 2 , as- 
suming a ballistic backward separation within a friction 
timescale. 

This picture for the rms relative velocity can be easily 
generalized to understand the behavior of the full PDF 
as a function of St. Consider a pair of particles at a dis- 
tance r at time 0, and suppose their relative velocity is 
w. Applying the above physical picture to this partic- 
ular particle pair, the relative speed, w, is estimated as 
w ~ Au(r p ), where Au is the flow velocity difference the 
two particles "saw" at r = — r p and r p is the primary 
distance of this particle pair. We have used w and Au to 
represent either the radial or tangential component. The 
generalized picture suggests that the particle relative ve- 
locity samples the PDF of the flow velocity difference in 
a certain way. An immediate implication of this picture 
is that the relative velocity would inherit intermittency 
of the turbulent flow. Assuming a ballistic separation 
again, r p is estimated by (r 2 + C,w 2 t 2 ) 1 / 2 , where £ ~ 3 
corresponds to the difference between the 3D separation 
velocity of the particle pair and the radial or tangen- 
tial component. We point out that, for particles with 
0.8 <, St <, 6.2, it may not be valid to assume the con- 
tribution to the particle relative speed is dominated by 
the ballistic separation phase. As discussed in §6.1, the 
Richardson phase may provide a crucial contribution for 
these particles (see Fig. [5]). However, using the ballis- 
tic assumption to estimate r p would be sufficient for a 
qualitative understanding of the relative velocity PDF. 

The above argument provides a satisfactory explana- 
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tion for our simulation results for the relative speed PDF, 
P(w, t p ), of particles with r p < Tl. At St -^i 1, the pri- 
mary distance r p (= (r 2 + C,w 2 T p ) x / 2 ) for particle pairs 
at the inner part of the PDF (i.e., w ~ 0) would be 
close to r. As a result, the central PDF follows the PDF, 
P u (Au, £), of the flow velocity difference at £ = r, as ob- 
served in the left panels of Figs. fT2l and [T3l At the tails 
of P(w,t p ), r p is larger, and w samples the flow veloc- 
ity PDF P u (Au,£) at larger £. This implies that higher 
tails broaden faster because P u (Au,£) is wider at larger 
£. The effect may be viewed as a "self-amplification" of 
the PDF tails. The tail amplification makes the overall 
shape of P(w,t p , St) at St & 1 considerably fatter than 
the PDF of tracers. As St increases, r p becomes larger at 
the same value of w, and the "amplification" proceeds to- 
ward the inner part of the PDF, as seen in the left panels 
of Figs. [T2l and [T3l The overall PDF broadening appears 
to be driven by the tail amplification. The amplification 
in the far PDF tails of St ^ 1 particles actually corre- 
sponds to the effects of slings or caustic formation. This 
is because the tail of P(w,t p ) is associated with local 
flow regions with large velocity gradients, which are in- 
deed where the slings or caustics are expected to occur. 
The tail amplification thus corresponds to the caustic 
contribution to the particle collision rate in the model 
of Wilkinson et al. (2006). On the other hand, the cen- 
tral part of the PDF represents the continuous (or S-T) 
contribution (Wilkinson et al. 2006). 

As St increases above 1, the range of the central PDF 
that follows P u (Au,£) becomes narrower, and the outer 
parts continue to get more extended. As discussed in 
§6.3.1, for St ~ 3.11, the PDF tails show slope changes 
as they extend to high values of \w\. This is because dif- 
ferent parts of the relative velocity PDF samples the flow 
PDF, P u (Au, £), at different length scales. As the fatness 
of P u (Au,£) decreases with increasing £ (see Appendix 
B and Fig. [HI), the shape of P{w, r p ) at higher tails be- 
comes thinner. This thinning trend occurs at smaller 
values of \w\ for particles with larger r p . The trend ex- 
plains why the overall shape of the PDF becomes less fat 
as St increases above ~ 1 . Note that the central cusp for 
3.11 & St <, 24.9 keeps a sharp shape, corresponding to 
P n {Au,£) at small £. 

The fact that the broadening of the PDF starts from 
the tail amplification is not captured by the PP10 model 
for the rms relative velocity in §3.2. The model only 
considers the 2nd-order moments of the flow velocity in- 
crement, the particle separation and the particle relative 
velocity, and thus implicitly assumes that the PDF shape 
does not change significantly with St, or the change oc- 
curs primarily at central part of the PDF. This gives rise 
to uncertainties in the prediction for the rms relative ve- 
locity because the PDF P(w, St) is found to be very fat 
especially for St ~ 1. Even the far tails give consider- 
able contribution to the variance. The tail amplification 
also provides evidence for a positive correlation between 
the fluctuations of the flow velocity increment and the 
particle separation. The tails of P(w,t p ) correspond to 
the PDF tails of both the flow velocity difference, Au, 
and the the primary distance, r p . In other words, in flow 
regions with Au larger than its rms value, the backward 
separation of particles is also faster than the rms separa- 
tion rate. As discussed in §3.2 and §6.1, the PP10 model 



neglects this correlation, and thus tends to underestimate 
the rms relative speed, particularly at intermediate St. 
The effect of the PDF tail amplification on the variance 
of the relative speed may be incorporated in the PP10 
model if the correlation between Am and r p are properly 
accounted for. 

In principle, if the PDF, P u (Au, £), of the flow velocity 
increment as a function of the length scale £ is provided, 
one can derive the PDF of the relative velocity as a func- 
tion of the particle size. For illustration, we consider a 
simple example, where the flow velocity is assumed to be 
exactly Gaussian. In that case, P n (Au,£) is written as, 



P u (Au,£) 



1 



exp 



Au 2 \ 
2W)J 
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where S(£) is the structure function or the variance of Au 
at I. To estimate the PDF, P(w, r p ), of the particle rel- 
ative speed w, we ask the question what the probability 
is for two nearby particles to see a flow velocity differ- 
ence of w at time r ~ — r p . The probability is roughly 
estimated by cx P u (w 7 r p ). Using eq. (|3"5"|) for P u and a 
ballistic particle separation for r p , we have, 



P(w, t p ) oc exp 



w 



2£S({r 2 + C,w 2 t 2 ) 1 / 2 ) 



(36) 



where it is assumed all the uncertainties in the rough 
estimate can be absorbed in a parameter £. If r p is in 
the inertial range of the flow, we may apply the Kol- 
morgorov scaling for S(£), i.e., S(£) oc £ 2 ^ 3 . Further as- 
suming that £ is independent of w, we find that eq. (|36[) 
corresponds to a stretched exponential with a = 4/3 (see 
eq. ([34]) ) at w r/r p . This suggests that the relative 
speed of inertial-range particles would be non-Gaussian 
even if the flow statistics were exactly Gaussian. This 
non-Gaussianity originates purely from the particle dy- 
namics, and is thus distinct from the non-Gaussianity in- 
herited from the intermittency of the turbulent flow. In 
other words, we identified two sources, namely, the flow 
intermittency and the particle dynamics, that contribute 
to the non-Gaussianity of the particle relative velocity. 

The predicted stretched exponential with a — 4/3 was 
found to satisfactorily fit the PDF tails of St = 24.9 par- 
ticles (see dashed lines in the left panels of Figs. 1121 and 
I13p . In this particular case, the two assumptions made 
in the derivation, i.e., Gaussianity and Kolmogorov scal- 
ing of the flow velocity, are both satisfied. However, we 
note that these assumptions are strong, and the validity 
of the 4/3 stretched exponential is quite limited. In fact, 
the prediction applies only to particles around the peak 
Stokes number, St m ~ 30. As discussed in §6.1, for parti- 
cles with St ~ Stm in our simulation, the typical primary 
distance is around 200?y. From Fig. [HI we see that, above 
this length scale, the PDF of the flow velocity increment 
is close to Gaussian. Therefore, the Gaussian assump- 
tion made in eq. (|36]) is valid for St ^ St m . Also, Fig. [3] 
shows that 2OO77 is toward the end of but still within the 
inertial range of the flow, meaning that, only for par- 
ticles with St <j St m , can one apply the Kolmorgorov 
scaling around the primary distance, r p . These suggest 
that the two assumptions are simultaneously met only 
at St ~ St m . Our finding that the 4/3 stretched expo- 
nential fits the PDF tails of St = 24.9 particles confirms 
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Fig. 14. — The normalized PDF of the radial relative velocity, w T , as a function of the Stokes number St and the particle distance, r. All 
PDFs are normalized to have unit variance. The normalized relative speed is defined as = w r /(w^ The solid, dashed and dotted 
lines correspond to particle distance r = 1,0.5 and 0.25?;, respectively. The left panel plots the PDFs for particles with St < 1.55, while 
the right panel shows the results for larger particles with St > 3.11. In each panel, the bottom lines (i.e., St = 0.1 and St = 795) show the 
actual PDF values, and, for clarity, the upper lines for each larger St are shifted upward by a factor of 16. 



the validity of our physical picture. But this stretched 
exponential form does not apply to the central part of 
the PDF, where both assumptions beak down. We find 
that the 4/3 stretched exponential can also acceptably fit 
the PDF tails for St = 12.4 particles, but not for other 
particles. 

We next consider large particles with St <; St m . The 
friction time of these particles is much larger than TL, 
and, accounting for the effect of the memory time of 
the flow velocity, the relative velocity of a given particle 
pair is roughly estimated by w ~ Au(r p )(Ti^/T p ) 1 / 2 (see 
§3.2.4). Due to the large friction time, r p is typically 
comparable to or even larger than the integral length 
scale, L, of the turbulent flow, meaning that P u (Au,£) 
at £ ~ r p is close to Gaussian. Since the flow velocity 
PDF "seen" by St St m particles is typically Gaussian, 
the PDF shape for their relative velocity appears to be 
simpler than particles with intermediate r p (see Fig. 1121 
an [T3")) . Using eq. (f3"5| and the same analysis that led to 
eq. ([361), we find P(w,t p ) oc exp(-(w 2 T p )/(2£S(r p )T h ) 
for St <; St m . As r p increases, r p increases, and S(r p ) 
becomes less dependent on r p or w. As a consequence, 
the shape of the relative velocity PDF becomes less fat. 
In the limit r p — > oo, S(r p ) would become a constant, 
2u' 2 , and P(w,t p ) finally approaches a Gaussian PDF 
with variance cx u' 2 Tl/t p . As observed in Figs. IT21 and 
[T3l a nearly Gaussian PDF is indeed observed for the 
largest particles in our simulation. 

6.3.3. The normalized PDF of the radial component 

To see the shape of the PDF more clearly, we show 
in Fig. [14] the PDF of the radial component normalized 
to have unit variance. The radial relative speed is nor- 
malized to its rms value, (w 2 ) 1 / 2 . The solid, dashed and 
dotted lines lines are normalized PDFs at r = lrj, 0.5 and 
0.25/7, respectively, and curves of different colors corre- 
spond to different St. The bottom curves in the left and 
right panels plot the actual PDF values for St = 0.1 and 



St = 795, respectively. For clarify, we shifted the PDF 
curves upward by a factor of 16 for each larger St in 
the left panel or each smaller St in the right panel. The 
asymmetry of the PDFs at St ^ 1.55 is clearly seen. 

As St increases from 0.1 to 1.55, the central part of 
the normalized PDF P(wr, St) becomes sharper. Note 
that, before normalization, the innermost part of the 
PDF is essentially the same for all particles in the range 
0.1 < St < 1.55 (see Fig. QJ]). Since the rms of w r in- 
creases with St, normalizing u> r by its rms tends to make 
the central part of the PDF sharper. This sharpening 
explains the decreases of the ratio (|w r |)/(w 2 ) 1 / 2 with 
increasing St for St < 1.55, as shown in the inset of the 
left panel in Fig. [KjJ At St > 3.11, the central cusp in 
the normalized PDF becomes so narrow that both the 
1st and 2nd order moments of the PDF are dominated 
by the outer parts of the PDF. Because the outer parts 
are less fat with increasing St, the ratio (|w r |)/(ix; 2 ) 1 / 2 
starts to increase at St ^ 3.11 (see Fig.^UJ. The behav- 
ior of (\w T \) I (if 2 ) 1 / 2 with St is consistent with the overall 
fatness of the PDF, i.e., it is smaller for fatter PDFs. 

It is interesting to note that, as St increases from 0.1 
to 0.39, the PDF, Piw^, St), at small to intermediate uf T 
in the right wing decreases. This corresponds to the de- 
crease of (w 2 ) 1 / 2 in the St range from 0.1 to 0.39 (see the 
left panel of Fig. [8j ■ The physical reason is that the right 
wing corresponds to separating particle pairs, and the 
particle distance decreases toward the near past. This 
leads to a decrease in the primary distance, r p , for sepa- 
rating pairs with small to intermediate w^, as St increases 
from 0.1 to 0.39. For the far right tail with large w T , the 
particle pairs may quickly move past each other, and 
their distance starts to increase within a friction time in 
the past. This explains why the far right tail of St = 0.39 
becomes slightly broader than that of St = 0.1 particles. 
As St increases to 0.78, the PDF at immediate to large 
Wx in the right wing is significantly amplified, while the 
effect of the "initial" decrease of the particle distance 
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is still visible at small positive uj. For St > 0.78, the 
particle memory is longer, and the "initial" separation 
phase does not cause a significant difference in the pri- 
mary distance r p for separating and approaching particle 
pairs. The two wings become symmetric at St ^> 1.55. 

For particles with St ^ 3.11, the normalized PDF has a 
dependence on r. As r decreases, the central part of the 
normalized PDF becomes sharper, and the outer parts 
become slightly broader, leading to a fatter PDF shape 
at smaller r. This can be understood as follows. Before 
the normalization, the central part of the PDF follows the 
flow velocity difference at r, and its width thus decreases 
linearly with r. On the other hand, the dependence of 
the tails on r is weaker because the primary distance, 
r p , for particle pairs in the outer parts of the PDF has a 
larger contribution from backward separation. Also, as 
r decreases, the contribution from the outer parts to the 
variance increases. Consequently, normalizing the PDF 
to unit variance gives a sharper central part and broader 
tails at smaller r. At St J> 6.21, the PDF is essentially 
independent of r. For these larger particles, r p is mainly 
contributed by the backward separation even for pairs 
lying close to the central part of the PDF, and thus the 
PDF is independent of r for w r in any range. The in- 
dependence of the normalized PDF shape for St <, 6.21 
explains the decrease of the ratio (\w T \) / (w 2 ) 1 / 2 with de- 
creasing r, as seen in the inset of the left panel of Fig. [10] 
(see §6.2). 

6.3.4. The tangential PDFs for approaching and separating 
particle pairs 

In Fig. [15l we compare the PDFs of the tangen- 
tial component of the relative velocity for approaching 
(P(w t \w T < 0,St)) and separating (P(u>t|u> r > 0, Si)) 
particle pairs. For St <, 6.2, the PDF of approaching par- 
ticles is broader than the separating ones, consistent with 

our earlier results for the rms relative speeds, (w 2 )^ 2 , 
shown in the left panel of Fig. [SJ Again, this is because, 
for a given "initial" value r, the distance of approach- 
ing particles was larger in the near past than the sep- 
arating ones. Therefore, the PDF of the relative veloc- 
ity for approaching pairs samples the PDF, P u (Au,£), 
of the flow difference at larger £. Since the width of 
P u (Au,£) increases with I, P(wt\w r < 0, St) is expected 
to be broader than P(w t \w T > 0, St). At St £ 12.4, the 
PDFs for approaching and separating pairs are almost 
equal. For these larger particles, the primary distance at 
r ~ — r p becomes insensitive to the initial conditions at 
r = 0. Although P(w t \w T > 0, St) for separating parti- 
cle pairs is not relevant for particle collisions, a compar- 
ison of P(wt\w T < 0, St) with P(wt\w T > 0, Si) provides 
an interesting illustration for the role of the backward 
separation of particle pairs in determining their relative 
velocity. 

6.3.5. The normalized PDF of the 3D amplitude 

Fig. [16] plots the PDF of the 3D amplitude, \w\, of the 
relative velocity for approaching particle pairs at r = lr/. 
For each St, the amplitude |to| is normalized to the rms 

value (w 2 ) 1 / 2 , i.e., |to| = \w\/ (w 2 )]/ 2 , so that all the 
PDFs shown here have unit variance. The rms of the 
3D amplitude, (w 2 )_ , for approaching pairs has been 
shown in Fig. [8] For a particle pair colliding with a rel- 
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ative velocity to, the total collision energy in the center- 
of-mass frame is imp|«;| 2 where is the reduced mass. 
Therefore, the PDF of the 3D amplitude can be easily 
converted to the PDF of the collisional energy. As men- 
tioned earlier, for particles of a given size, the PDF of 
| to | determines the fractions of collisions that result in 
sticking, bouncing or fragmentation. 

The left and right panels of Fig. [T5] show simulation 
results for small (St <, 1.55) and large (Si 3.11) parti- 
cles, respectively. The thin dashed line in the left panel 
corresponds to approaching tracer particles (Si = 0), 
while the dashed line in right panel is the normalized 
PDF of a Gaussian vector with three independent com- 
ponents of equal variance. The PDF for tracer particles 
in the left panel is already highly non-Gaussian, as can 
be seen from a comparison with the dashed line in the 
right panel. In the left panel, we see that the degree of 
non-Gaussianity increases as Si increases from to ~ 1. 

At larger Si, the PDF peaks at smaller |to|, indicating a 
larger probability of collisions with lower collision energy. 

The PDF around the rms value (i.e., |uj| ~ 1) decreases 
with increasing Si, and more probability is distributed 

toward smaller and larger values of \w\. This corresponds 
to the sharpening of P(w r ,St) and P(w t \w r < 0,St) in 
the central part and the broadening of the tail parts in 
this Si range (see Figs. [121 [TBI and Ej). 

The trend is reversed as Si increases further above 
St ~ 3.11. The peak of the PDF moves back to around 

the rms value, \w\ ~ 1, at St ^ 49.7. The PDF eventu- 
ally approaches Gaussian in the limit r p 3> Tl- However, 
note that, even for Si = 795 particles, the PDF shows 
a significant difference from the Gaussian distribution at 
small collision speeds. The strong non-Gaussianity in the 
amplitude of the collision velocity has interesting impli- 
cations for the growth and evolution of dust particles in 
planetary disks. The shape of the PDFs shown in Fig. 
[16] suggests that there are more collisions with extremely 
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large (\w\ ^> 1) or small (\w\ ^S> 1) relative speeds than 
estimated from a Gaussian distribution. This may in- 
duce interesting and complicated effects to the model- 
ing of the collision process. The high probability for 
collisions with small relative speed would favor sticking, 
while there are also more collisions that would result in 
fragmentations. The competition of the two opposite ef- 
fects would determine whether the non-Gaussian PDF 
of the collision velocity accelerates or slows down the 
particle growth. A coagulation model incorporating the 
non-Gaussian statistics of the collision speed would give 
a more realistic prediction for the evolution of the size 
distribution of dust particles in protoplanetary disks. 

We find that, for Si in the range from 0.78 to 6.22, the 
PDF shows an extended power-law range at intermediate 

values of \w\. For example, at St — 0.78, the PDF goes 

like |mj| -2 - 4 in the range 0.5 < \w\ < 4. The slope of 
the PDF in the power-law range becomes shallower with 
increasing St. For St = 1.55, St = 3.11 and St = 6.22, 
the power-law exponent of the PDF in the intermediate 
range is —1.8, —1.3, and —0.8, respectively. 

The PDF of the amplitude can be easily computed 
from the PDFs of the radial and tangential components, 
if the three components are completely independent. In 
that case, one may obtain fitting functions for P(|u;|) 
using the fitting functions discussed earlier for the PDFs 
P(w T ,St) and P(wt\w T < 0, Si), of colliding particle 
pairs. Alternatively, one could directly fit the PDF of \w\ 
with simple function forms or tabulate it as a function of 
St. We defer this to a later paper using simulations at a 
higher resolution. 

7. CONCLUSIONS 

We investigated the turbulence-induced relative veloc- 
ity of inertial particles using both theoretical modeling 
and numerical simulations. We conducted a 512 3 simu- 
lation of a weakly compressible turbulent flow with an 
rms Mach number of 0.1, and evolved 14 species of in- 
ertial particles with friction timescales, r p , covering the 



entire scale range of the flow. The Stokes number, Si, of 
the particles spans about 4 order of magnitudes, ranging 
from 0.1 to 800. We used the simulation to test the theo- 
retical model for the rms velocity of inertial particles by 
Pan & Padoan (2010), in the case of identical particles 
(equal friction times). We also examined the particle 
collision kernel and the probability distribution of the 
relative velocity as a function of St. In the following, we 
summarize the main results of this work. 

1. We introduced the formulation for the particle rel- 
ative velocity by Pan & Padoan (2010), which re- 
veals an insightful physical picture. The relative 
velocity of two nearby identical particles is deter- 
mined by the memory of the flow velocity difference 
along their trajectories in the past, and hence de- 
pends on the separation behavior of particle pairs 
backward in time. We adopted a two-phase sep- 
aration behavior consisting of a ballistic and a 
Richardson phase, and showed that the model pre- 
dicts a Si 1 / 2 scaling for inertial-range particles in 
turbulent flows with an extended inertial range. 
The model can correctly reproduce the expected 
behaviors in the extreme limits of large and small 
particles. We found that the model prediction for 
the rms relative velocity is in good agreement with 
the simulation results. The physical picture also 
provides a successful explanation for the behavior 
of the probability distribution function of the rela- 
tive velocity as a function of St. 

2. To improve the understanding of the inertial par- 
ticle statistics, we analyzed both Lagrangian and 
Eulerian temporal correlation functions in the sim- 
ulated flow. While the flow velocity along the tra- 
jectory of a small particle is close to Lagrangian, 
the velocity seen by a large particle may be better 
approximated by Eulerian. The Eulerian and La- 
grangian correlation timescales, Te and X"l, were 
found to be close to each other, with Te slightly 
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larger (by 10%). Our model predictions for both 
1-particle rms velocity and the rms relative speed 
between two identical particles depend mainly on 
the correlation timescale and are insensitive to the 
function form of the temporal correlation. This 
provides a validation for using the Lagrangian cor- 
relation function form for all particles. 

3. Our simulation shows that, in the small particle 
limit (St <C 1), the 3D rms relative velocity of par- 
ticle pairs at a distance r = rj is constant, ~ u n / >/3, 
consistent with the Saffman- Turner prediction. It 
starts to rise at St ^ 1, and peaks for particles with 
r p ~ 2Tl, corresponding to St = St m ~ 30 in our 
simulated flow. As expected, the relative velocity 
scales with St as St -1 / 2 in the limit r p 3> Tl (or 
St » Stm)- The PP10 model with reasonable pa- 
rameters provides an excellent fit to the simulation 
data for the 3D rms. The maximum relative speed 
at St m is twice smaller than the rms velocity of the 
turbulent flow, indicating a factor of ~ 2 overesti- 
mate by the commonly-used models of Volk et al. 
(1980) and its later developments (e.g., Markiewicz 
et al. 1991; Cuzzi & Hogan 2003; Ormel & Cuzzi 
2007). 

The rms relative speed of small particles (St ^ 6) 
has a dependence on r. The dependence for the 
smallest particles (St — 0.1 — 0.2) in our simulated 
flow at r < rj was found to be slower than the 
linear scaling predicted by the Saffman- Turner for- 
mula, suggesting considerable contributions from 
the sling events or caustic formation. In the con- 
text of the PP10 model, these events provide a con- 
tribution to the backward particle separation. For 
larger particles with St <; 6, the backward separa- 
tion at a friction timescale ago becomes insensitive 
to the initial distance r, and the rms relative speed 
is thus independent of r. 

The rms relative speeds in the radial and tangential 
directions are nearly equal for all St ^ 0.1 particles 
at r ^ \rj. For St -C 1 particles, this is in contrast 
to the Saffman- Turner formula, which predicts the 
tangential rms is larger than the radial rms by 
v/2. This near equality suggests the direction of 
the relative velocity is almost completely random- 
ized with respect to the particle separation r. The 
randomization is caused by the derivation of parti- 
cle trajectories from the fluid elements and/or the 
stochastic backward separation of particle pairs. 

For St 6 particles, approaching pairs that may 
lead to collisions were found to have a larger rela- 
tive speed than separating ones. The physical rea- 
son is that the distance between approaching parti- 
cles is larger than separating ones in the near past. 
This implies that the collision speed of St ^ 6 par- 
ticles is larger than the overall relative velocity es- 
timated from all particle pairs at a given distance. 

4. We computed the particle collision kernel from the 
simulation data using both spherical and cylindri- 
cal formulations. In particular, we accounted for 
the effect of turbulent clustering. It was found that 
the two formulations give nearly equal predictions 



for all particles. The collision kernel per unit cross 
section shows an abrupt rise as St increases toward 
1, which may be viewed as an activation process 
corresponding to the formation of caustics. As St 
increases from ~ 1 to Stm, the normalized collision 
kernel is roughly constant, increasing only slightly 
by 50%. It finally decreases as St^ 1 / 2 for large 
particles with t p 3> Tl. The normalized kernel is 
independent of the particle distance, r, for St <; 1, 
due to the cancellation of the scalings of the ra- 
dial distribution function and the collision velocity 
with r. On the other hand, for St <, 1 particles, 
the kernel shows a dependence on r, which needs 
to be accurately evaluated for applications to pro- 
toplanetary disks, where the dust grains are much 
smaller than the Kolmorgorov scale. At sufficiently 
small r, the effect of caustics may provide a dom- 
inant and r— independent contribution. Existing 
coagulation models for dust particle growth do not 
capture these important features. We will provide 
fitting functions for the collision kernel as a func- 
tion of St in a separate study. 

5. The probability distribution function for the par- 
ticle relative velocity is highly non-Gaussian, ex- 
hibiting extremely fat tails. For small particles 
with St <, 1, the effects of the particle memory and 
the backward separation lead to a self-amplification 
starting at the far tails, which corresponds to slings 
or caustic formation. As St increases, the amplifi- 
cation becomes stronger and proceeds toward the 
inner part of the PDF, causing an increase in the 
fatness of the overall PDF shape. On the other 
hand, as St increases above ~ 1, the PDF shape 
becomes less fat. For these larger particles, the 
relative velocity samples the PDF, P u (Au,£), of 
the flow velocity increment, Au, at larger scales, 
I. Since the fatness of P u decreases with increasing 
£, the PDF of the particle relative velocity keeps 
thinning with increasing St. At a particle distance 
of r = lry, the PDF shape is fattest at St ^ 1, with 
a kurtosis of ~ 30. For St <J 6, the PDF shape also 
depends on r and becomes even fatter at r & lr/. 

Using the physical picture of the PP10 model, we 
identified two sources of non-Gaussianity for the 
particle relative velocity: the inheritance of in- 
termittency from the turbulent flow and an in- 
trinsic contribution from the particle dynamics. 
We also predicted a 4/3 stretch exponential PDF, 
cx cxp(— (|i/j|//3) 4 / 3 ), for inertial-range particles in 
an exactly Gaussian velocity field with Kolmor- 
gorov scaling. This 4/3 stretched exponential is ac- 
tually observed in the PDF tails for particles with 
St ~ Stm (or t p ~ 2TL), confirming the validity of 
our physical picture. 

The PDF of the particle collision velocity is ex- 
pected to play a crucial role in the growth of dust 
particles in protoplanetary turbulence, as it deter- 
mines the fractions of collisions resulting in stick- 
ing, bouncing and fragmentation. The PDF of the 
3D amplitude of the relative velocity shows much 
higher probabilities of extremely small and large 
collision speeds than estimated from a Gaussian 
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PDF. Only for particles with r p <; 50Xl does the 
PDF approach Gaussian. The highly non-Gaussian 
nature of the collision velocity needs to be incorpo- 
rated into dust coagulation models for protoplane- 
tary disks. 



In this work, we have focused on the monodisperse case 
with identical particles. A systematic analysis for the re- 
late velocity of different particles will be conducted in a 
follow-up paper. Future simulations at higher resolutions 
will further improve our understanding of the problem. 



For example, a 1024 3 simulation would help verify the 
existence of the predicted St 1 / 2 scaling for inertial-range 
particles by various models, which is not yet confirmed 
numerically. A computationally more demanding sim- 
ulation with a much larger number of inertial particles 
would allow us to examine the particle statistics at small 
scales, r -C r], and to quantify the r-dependence of the 
collision kernel and the caustic contribution for St <, 1 
particles. The results of these future studies will signifi- 
cantly improve the formulation of coagulation models to 
compute the evolution of dust particles in protoplanetary 
disks. 



APPENDIX 



A SEPARATION OF TRACER PARTICLE PAIRS 

Our model for the relative velocity of inertial particles depends on the particle pair dispersion backward in time 
(§3.2.3). We adopted a two-phase behavior, consisting of a ballistic and a Richardson phase. To constrain the 
Richardson constant, g, in the latter phase, we study the dispersion of tracer particles in our simulated flow, and take 
the measured value of g as a reference for inertial particles. 

In Fig. 1X3 the three solid lines from bottom to top show the backward-in-time separation of tracer particle pairs 
with "initial" distance r = I77, 2rj and 4?y, respectively. We subtracted r 2 from the separation variance, d 2 {r). The 
"initial" time is set to 0, and r is negative for the backward separation. As seen in Fig. 1171 the particle separation at 
small |t| shows a ballistic behavior with d 2 (r) — r 2 increasing as r 2 . The physical origin of the ballistic behavior is 
that the velocity at which two tracer particles separate is determined by the flow velocity difference across the particle 
distance and thus remains roughly constant before the particle distance becomes significantly larger than the initial 
value. This ballistic phase of the tracer particles is physically different from that of inertial particles discussed in 
§3.2.3. For inertial particles, the duration and the separation speed of the ballistic phase depend on the particles' 
memory timescale. But for tracer particles, the ballistic phase is determined purely by the initial particle distance. 
For r < 4?y, the ballistic phase lasts for a few Kolmorgorov timescales. 

The Richardson separation behavior (the |t| 3 scaling) is observed at large |r| after the particle separation enters the 
inertial range of the flow. In the bottom solid curve for r = I77, a |r| 3 scaling exists in a very limited range. A rough 
estimate of the Richardson constant in that range gives g ~ 0.5. Consistent with previous studies (e.g., Sawford et al. 
2008), the time range that exhibits the Richardson scaling becomes broader as r increases to 4rj. This allows a more 
accurate estimate of g, and we find g ~ 1.2 for r = 4?/, consistent with the experimental results of Berg et al. (2006). 
Similar to Sawford et al. (2008), the measured g has a dependence on r. If the inertial range of the flow is considerably 
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broader and the Richardson separation behavior exists in a larger time range, then the curves for different initial 
distances are expected to converge at sufficiently large time lags, with g eventually approaching a universal value. 

The dashed line in Fig. [17] plots the forward-in-time pair dispersion of tracer particles at r — Arj. A comparison 
with the top solid line shows that the forward separation is slower than the backward separation. For r — 4?y, g is 
estimated to be 0.5 in the forward separation, about twice smaller than the value (1.2) for the corresponding backward 
separation. This is consistent with the result of Berg et al. (2006). A similar trend is also found for r = \r\ and 1r\. 
The faster backward separation will be explained in Appendix B based on the PDF of the longitudinal difference, Ait r , 
of the flow velocity. In our model for the relative velocity of incrtial particles, it is the backward separation that is 
relevant, and the purpose of showing an example for the forward separation in Fig. [IT] is to illustrate the difference 
between the forward and backward separations. 

The Richardson constant, g, in the tracer-like phase of inertial particle separation may be different from tracers. 
However, it is reasonable to assume that g for the backward separation of inertial particles lies in a similar range. Like 
tracers, the value of g for inertial particles may also depend on the initial distance, r. In §6.1, we adjusted value of g 
in our model to obtain best fits to the simulation data for the particle relative velocity at different r. 

B THE PDFS OF THE TURBULENT VELOCITY FIELD 

In this Appendix, we analyze the probability distribution functions (PDF) of the flow velocity increments at different 
length scales. As discussed in §6.3, the statistics of the flow velocity increment leaves a signature in the PDF of the 
relative velocity of inertial particles. We measured the PDFs, P u (Au T ,£) and P u (Au t ,£), of the longitudinal and 
transverse increments as a function of the length scale, £. Similar to the computation of the structure functions 
in §4.2, the PDF measurement also uses the velocity differences along the three directions of the simulation grid. 
The variances of P u (Au T ,£) and P u (Aut,£) correspond to the longitudinal (S\\(£)) and transverse (S nn (£)) structure 
functions, which have been shown in Fig. [3] (see §4.2). Clearly, the PDFs are wider at larger scales as the structure 
functions increase with £. Furthermore, P u (Au t ,£) is wider than P u (Au T ,£) because S nn (£) > Su(£) at all £ (see Fig. 

El. 

To better see the PDF shape as a function of £, we normalized the PDFs at each scale to have unit variance. 

' 1/2 

The radial and transverse velocity increments are normalized to their rms values, i.e., Au r (£) = Au T {l)/S^ (I) and 

Au t {£) = Aw t W/5nn 2 (^). The left panel of Fig. [TBI shows the normalized PDF, P u (Au r ,£), of the radial increment. 
Except at the largest scales, the PDF is negatively skewed. For inertial-range scales, this can be understood from 
Kolmorgorov's 4/5 law, (Au r (£) 3 ) = — | et, which indicates a negative skewness for the PDF of Ait r . The connection of 
the 3rd order moment of Au r to the energy dissipation rate suggests that the skewness originates from the dissipative 
nature of turbulence. The skewness of the PDF of Aw r also provides an explanation for the faster backward separation 
found in Appendix A. The left and right tails of P u (Au r , £) correspond to tracer pairs receding from each other backward 
and forward in time, respectively. The broader left tail of the PDF thus suggests that the backward separation of 
tracer particles is faster than the forward case. The argument also implies that the time-asymmetry (or irreversibility) 

of the tracer pair separation is related to the dissipative nature of a turbulence system. Unlike P u (Au r ,£), the PDF, 




I! 




Fig. 18. — The normalized PDFs of the flow velocity increment in the radial (left panel) and transverse (right panel) directions as a 
function of length scale, i. The normalized velocity increments are defined as Au r (I) = Aur(£)/sl/ 2 (£) and Au t (£) = Au t (£) / 'S^l? '(£) . The 
top line in each panel plots the exact PDF at £ = 1.7?) (the cell size), and, for clarity, the PDF at each larger £ is shifted downward by a 
factor of 4. Except at the largest scales, the PDF, P u (Au T , £), of the radial increment has a negative skewness, whereas P u (AMt,^) for the 
transverse increment is symmetric at all scales. The PDF tails are highly non-Gaussian or intermittent at small scales. With increasing £, 
the PDFs become less fat and approach Gaussian at the largest scales. In the right panel, the dashed and dotted lines for £ = 1.7rj are the 
normalized PDFs of the transverse increment conditioned on A« r < and A« r > 0, respectively. 
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P u (Au t , £), of the transverse increment in the right panel is symmetric about the origin, Au t = 0, at all I, as expected 
from statistical isotropy. 

Both P u (Au r , £) and P u (Au t , £ ) are close to Gaussian at the largest scales, € = 2II77 (1/4 box size) and 422?? (1/2 box 
size), of the simulated flow. This is consistent with the Gaussian 1-point statistics in fully developed turbulence. At 
smaller £, the PDFs become non-Gaussian, and the tails keep getting fatter with decreasing £, a phenomenon known 
as intermittency in turbulence theory (Frisch 1995). As mentioned in the text, we use the word "fat" (or "thin") 
specifically for the shape of the PDF, while "broad" (or "narrow") refers to the extension or width of the PDF. The 
smallest scale, 1.7ry, in the figure corresponds to the size, Ax, of the computational cell. The shape of the normalized 
PDF is expected to remain unchanged once £ becomes smaller than ~ 77. Physically, the viscosity acts to smooth the 
velocity field and makes it diffcrcntiable in the dissipation range, and consequently the velocity increment across any 
scale £ <, r\ is proportional to the local velocity gradient, whose PDF is fixed. In our simulation, the velocity field 
inside a computation cell is obtained by interpolation, and thus the PDF of the velocity difference below the cell size 
is controlled by the velocity gradient PDF at Ax. In §6.3, we showed that the trend of the PDF shape of the flow 
velocity difference with £ has interesting effects on the PDF of the relative velocity of inertial particles as a function 
of the particle inertia. 

The tails of P u (Au t ,^) for the transverse increment can be approximately described by stretched exponentials, P se 

(see eq. (J34j) in §6.3). At largest scales, P u (AM t ,^) are nearly Gaussian, and a = 2. With decreasing £, a decreases, 

corresponding to fatter tails. For example, the best- fit a for P u (Au t ) at £ = 26ry is ~ 1, and it further decreases to 

~ 0.72 at £ — I.777. Due to the asymmetry of the radial PDF, P u (Au r , £), one needs to obtain the fits separately for the 

left and right wings. Comparing the left wing of P(Au r , £) with P u (Aut,£), we see that their shape has a similar level 

of fatness at the same scale £. In fact, the best-fit a for the left tail of P u (Au r ,£) is very close to that for P u (Au t ,£). 

We also find that the best-fit values of a for the left and right tails of P u (Au T ,£) are close to each other, indicating 
that the two tails have a similar shape and differ only in the amplitude of fluctuations. 
To quantify the fluctuation amplitudes in the left and right wings of P u (Au r ,£), we define the 

variances in the two wings as (Aw 2 )_ = J_ oQ Au^P u (Au I ,£)dAu I /J_ Qc P u (Au I ,£)dAu T and (Au^.) + = 
J °° Av%P u {Au T ,t)dAu T / J oo P u (Au r ,£)dAu r . The definition of (Au^) T is similar to ( w T ) =p for the relative veloc- 
ity of inertial particles (§6.1.3.) We find that the ratio of (Au 2 )_ to (Au 2 ) + is ~ 1.47 at £ — 1.7t). This ratio decreases 
with increasing £, and reaches unity at the largest scales. The variances of the left and right wings of P u (Au T ,£) was 
used in the discussion on the relative velocity of approaching and separating particle pairs in the St <C 1 limit (see 
§6.1.3). We also considered the PDF of Au t conditioned on the sign of Ait r . We denote two conditional PDFs as 
P u (Awt|Aii r ^5 0,^) and their variances as (Au t 2 ) T = Au t 2 P u (Au t |Au r § 0,£)dAu t . At £ = 1.7/7, ( Au tV is found 
to be larger than (Ait 2 ) + by 28%. In the right panel of Fig. [TH1 the dashed and dotted lines show the normalized 
conditional PDFs, P u (A-u t |A-u r < 0,£) and P u (Aw t |Aw r > 0,f), at I = l.7rj. The shapes of the two conditional PDFs 

are close to that of the overall PDF, P u (Au t ,£). Again, the conditional variances and PDFs of Au t are useful to 
understand the differences in the relative velocity of approaching and separating inertial particle pairs with St <C 1 
(§6.1.3 and 6.3.4). 

REFERENCES 



Abrahamson, J. 1975, Chem. Eng. Sci. 30, 1371 
Ayala, O., Rosa, B. & Wang, L.-P. 2008, New J. Physics, 10, 
075016 

Balkovsky, E., Falkovich, G., & Fouxon, A. 2001, Phys. Rev. 
Lett., 86, 2790 

Bee, J., Biferale, L., Cencini, M., Lanotte, A. S., & Toschi, F. 

2009, arXiv: 0905.1192vl [physics.fly-dyn] 

Bee, J., Biferale, L., Lanotte, A. S., Scagliarini, A., & Toschi, F. 

2010, J. Fluid Mech., 645, 497 

Berg, J., Luthi, B., Mann, J. & Ott, S. 2006, Phys. Rev. E., 74, 
016304 

Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, 
11 

Blum, J. & Wurm, G. 2008, ARAA, 46, 21 
Brandenburg, A. & Dobler, W. 2002, Computer Physics 

Communications, 147, 471 
Chambers, J. E. 2010, Icarus, 208, 505 
Chiang, E. 2008, ApJ, 675, 1549 
Cuzzi, J. N. & Hogan R. C. 2003, Icarus,164, 127 
Cuzzi, J. N., Hogan, R. C, & Bottke, W. F. 2010, Icarus, 208, 518 
Cuzzi, J. N., Hogan, R. C, Paque, J. M., & Dobrovolskis, A. R. 

2001, ApJ, 546, 496 
Cuzzi, J. N., Hogan, R. C. & Shariff, K. 2008, ApJ, 687, 1432 
de Jong J., Salazar, J. P. L. C, Woodward, S. H., Collins, L. R., 

& Meng, H. 2010, Int. J. Multiphase Flow, 36, 324 
Dullemond, C. P. & Dominik, C. 2005, A&A, 434, 971 



Falkovich, G., Fouxon, A., & Stepanov, M. G. 2002, Nature, 419, 
151 

Falkovich, G., & Pumir, A. 2004, Physics of Fluids, 16, L47 
Falkovich, G. & Pumir, A. 2007, J. Atmos. Sci. 64, 4497 
Freytag, B., Allard, F., Ludwig, H.-G., Homeier, D., & Steffen, 

M. 2010, A&A, 513, A19 
Frisch, U. 1995, Turbulence. The Legacy of AN Kolmogorov 

(Cambridge University Press, Cambridge) 
Girimaji, S. S. & Pope, S. B. 1990, J. Fluid. Mech., 220, 427 
Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 
Gustavsson, K. & Mehlig, B. 2011, Phys. Rev. E, 84, 045304(R) 
Gustavsson, K., Mehlig, B., Wilkinson, M. & Uski, V. 2008, Phys. 

Rev. Lett., 101, 174503 
Gustavsson, K., Meneguz, E., Reeks, M. & Mehlig, B. 2012, New 

J. Physics, 14, 115017 
Guttler, C, Blum, J., Zsom, A., Ormel, C. W.,& Dullemond, C. 

P. 2010, A&A, 513, A56 
Helling, Ch., Jardine, M. & Mokler, F. 2011, ApJ, 737, 38 
Hogan, R. C, & Cuzzi, J. N. 2001, Phys. Fluids, 13, 2938 
Hubbard, A. 2012, 426, 784 

Ishihara, T., Gotoh, T. & Kaneda, Y. 2009, Annu. Rev. Fluid 
Mech., 41, 165 

Johansen, A., Andersen, A. C, & Brandenburg, A. 2004, A&A, 
417, 361 

Johansen, A., Klahr, H., & Henning, Th. 2011, A&A, 529, A62 



35 



Johanscn, A., Oishi, J. S., Low, M.-M. M., et al. 2007, Nature, 
448, 1022 

Johansen, A. & Youdin, A. N. 2007, ApJ, 662, 627 

Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75 

Kolmorgorov, A. N. 1962, J. Fluid. Mech., 13, 82 

Kruis, F. E. & Kusters, K. A. 1997, Chem. Eng. Comm., 158, 201 

Lee, A. T., Chiang, E., Asay-Davis, X., & Barranco, J. 2010, 

ApJ, 725, 1938 
Lundgren, T. S. 1981, J. Fluid Mech. Ill, 27 
Markicwicz, W. J., Mizuno, H. & Voelk, H. J. 1991, A&A, 242, 

286 

Maxey, M. R. 1987, J. Fluid Mech., 174, 441 

Monin, A. S. & Yaglom, A. M. 1975, Statistical Fluid Mechanics: 

Mechanics of Turbulence, vol. 2. MIT press. 
Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 
Ormel, C. W., Paszun, D., Dominik, C, & Ticlcns, A. G. G. M. 

2009, A&A, 502, 845 
Padoan, P., Jimenez, R., Nordlund, A. & Boldyrev, S., 2004, 

Phys. Rev. Lett. 92, 191102 
Pan, L. & Padoan, P. 2010, J. of Fluid Mech., 661, 73 (PP10) 
Pan, L., Padoan, P.. Scalo, J., Kritsuk, A. G., & Norman, M. L. 

2011, ApJ, 740, 21 
Pan, L. & Scannapieco, E. 2011, Phys. Rev. E, 2011, 83, 

045302(R) 

Pinsky, M. B. & Khain, A. P. 1997, J. Aerosol Sci. 28, 1177 
Pruppacher, H. R., & Klett, J. D. 1997, Microphysics of Clouds 

and Precipitation (Dordrecht: Kluwer) 
Rossow, W.B., 1978, Icarus, 36, 1 

Saffman, P. G. & Turner, J. S.1956 J. Fluid Mech., 1, 16 
Sawford, B. L. 1991, Phys. Fluids, 3, 1577 

Sawford, Yeung& Hackl]saw08 Sawford, B. L., Yeung, P. K. & 

Hackl, J. F. 2008, Phys. Fluids, 20, 065111 
Shaw, R. A. 2003, Annu. Rev. Fluid Mech., 35, 183 



Squires, K. D., & Eaton, J. K. 1991, Phys. Fluids, 3, 1169 
Sundaram, S., & Collins, L. R. 1997, J. of Fluid Mech., 335, 75 
Volk, H. J., Jones, F. C, Mornll, G. E.& Roeser, S. 1980, A&A, 
85, 316-325 

Wang, L.-P., Wexlcr, A. S., & Zhou, Y. 1998, Phys. Fluids, 10, 
2467 

Wang, L.-P., Wexlcr, A. S., & Zhou, Y. 2000, J. of Fluid Mech., 
415, 117 

Weidenschilling, S. J. 1980, Icarus, 44, 172 
Wilkinson, M. & Mehlig, B. 2005, Europhys. Lett., 71, 186 
Wilkinson, M. Mehlig, B., & Bezuglyy, V. 2006, Phys. Rev. Lett., 
97, 048501 

Williams, J. J. E. & Crane, R. I. 1983, Int. J. Multiphase Flow, 9, 
421 

Yeung, P. K., Pope, S. B. & Sawford, B. L. 2006, J. Turbulence, 
7, 58 

Youdin, A. N. 2011, ApJ, 742, 38 
Yuu, S. 1984, AIChE J., 30, 802 

Zaichik, L. I., & Alipchenkov, V. M. 2003, Phys. Fluids, 15, 1776 
Zaichik, L. I. & Alipchenkov, V. M. 2009, New J. Phys., 11, 
103018 

Zaichik, L. I., Simonin, O. & Alipchenkov, V. M. 2003, Phys. 

Fluids, 15, 2995 
Zaichik, L. I., Simonin, O. & Alipchenkov, V. M. 2006, Phys. 

Fluids,18, 035110 
Zhou, Y., Wexler, A. S., & Wang, L.-P. 2001, J. Fluid Mech., 433, 

77 

Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, 
A&A, 534, A73 

Zsom, A., Ormel, C. W., Guttler, C, Blum, J., & Dullemond, 
C. P. 2010, A&A, 513, A57 



